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Abstract 

The  continued  proliferation  of  affordable  RF  communication  devices  has  greatly 
increased  wireless  user  exposure  and  the  need  for  improved  security  to  protect  against 
spoofing.  This  work  addresses  various  Open  Systems  Interconnection  (OSI)  Physical 
(PHY)  layer  mechanisms  to  extract  and  exploit  RF  waveform  features  (“fingerprints”) 
that  are  inherently  unique  to  specific  devices  and  that  may  be  used  for  reliable  de¬ 
vice  classification  to  provide  hardware  specific  identification  (manufacturer,  model, 
and/or  serial  number).  Automatically  detecting,  identifying  and  locating  RF  commu¬ 
nication  devices  remains  a  challenging  technical  problem  and  consists  of:  1)  the  selec¬ 
tion  and  generation  of  fundamental  signal  characteristics  (amplitude,  phase,  and/or 
frequency),  2)  the  feasibility  and  repeatability  of  detecting  and  locating  the  start 
of  a  burst  using  selected  waveform  feature(s)  amidst  channel  noise,  3)  the  identi¬ 
fication  and  robust  extraction  of  distinguishable  hngerprints-features  that  uniquely 
characterize  the  unintentional  modulation  of  a  device,  and  4)  the  performance  of  sig¬ 
nal  classification  under  varying  channel  conditions  and  Signal-to-Noise  Ratio  (SNR). 
This  challenge  is  addressed  by  applying  a  Dual-Tree  Complex  Wavelet  Transform 
(DT-CWT)  to  improve  burst  detection  and  RF  fingerprint  classification. 

Two  burst  detection  techniques  are  analyzed  under  varying  channel  SNR  con¬ 
ditions,  the  Fractal  Bayesian  Step  Change  Detector  (Fractal-BSCD)  and  Traditional 
Variance  Trajectory  (VT).  Performance  of  both  techniques  are  consistent  with  per¬ 
fect  burst  location  at  higher  SNRs  (10  <  SNR  <  30  dB)  but  diverged  at  lower  SNRs 
(—3  <  SNR  <  10  dB).  Traditional  VT  performance  is  most  consistent  with  perfect 
results  for  6  <  SNR  <  30  dB,  under  performs  perfect  results  for  —3  <  SNR  <  6  dB, 
and  outperforms  Fractal-BSCD  considerably  for  —3  <  SNR  <  18  dB.  A  “De- 
noised  VT”  technique  is  introduced  to  improve  performance  at  lower  SNRs,  with 
denoising  implemented  using  a  DT-CWT  decomposition  prior  to  Traditional  VT  pro- 
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cessing.  This  proves  to  be  effective  and  provides  more  robust  burst  detection  for 
-3  <  SNR  <  10  dB. 

Performance  of  a  newly  developed  Wavelet  Domain  (WD)  fingerprinting  tech¬ 
nique  is  presented  using  statistical  WD  fingerprints  with  Multiple  Discriminant  Anal¬ 
ysis/Maximum  Likelihood  (MDA/ML)  classification.  The  statistical  fingerprint  fea¬ 
tures  are  extracted  from  coefficients  of  a  DT-CWT  decomposition.  Relative  to  pre¬ 
vious  Time  Domain  (TD)  results,  the  enhanced  WD  statistical  features  provide  im¬ 
proved  device  classification  performance.  Improvement  is  characterized  using  a  “gain” 
metric  defined  as  the  difference  in  required  SNR  in  dB  ( SNRWd  ~  SNRTD)  for 
the  two  techniques  to  achieve  a  given  classification  performance.  Accounting  for 
all  intra-manufacturer  and  inter- manufacturer  device  discrimination  scenarios,  the 
WD  technique  provides  2-7  dB  of  gain  for  80%  correct  classification  performance  at 
2  dB  <  SNRWd  <  11  dB.  Additional  performance  sensitivity  results  are  presented 
to  demonstrate  WD  fingerprinting  robustness  for  variation  in  burst  location  error, 
MDA/ML  training  and  classification  SNRs,  and  MDA/ML  training  and  classification 
signal  types.  For  all  cases  considered,  the  WD  technique  proves  to  be  more  robust 
and  exhibited  less  sensitivity  when  compared  with  the  TD  technique. 
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Application  of  Dual-Tree  Complex  Wavelet  Transforms 
to  Burst  Detection  and  RF  Fingerprint  Classification 


I.  Introduction 

This  chapter  introduces  the  dissertation  research  and  its  documentation.  The  moti¬ 
vation  for  conducting  the  research  is  first  provided  in  Section  11.11  which  includes  the 
Operational  Motivation  factors  in  Section  11.1.11  and  Technical  Motivation  factors  in 
Section fl.  1.21  This  is  followed  by  a  summary  of  Research  Contributions  in  Section [L2l 
which  provides  a  relational  mapping  in  Table  11.11  to  highlight  contributions  of  this 
work  relative  to  what  had  been  previously  accomplished.  The  chapter  concludes  with 
a  Dissertation  Overview  in  Section  11.31 

1.1  Research  Motivation 

1.1.1  Operational  Motivation.  The  continued  proliferation  of  inexpensive 
wireless  Radio  Frequency  (RF)  devices  provides  worldwide  communication  connec¬ 
tivity  to  virtually  every  individual.  Within  a  geographically  localized  region,  the 
fundamental  emissions  from  these  devices,  i.e.,  the  intentionally  radiated  emissions 
designed  to  support  the  intended  purpose,  may  be  remotely  intercepted  by  unintended 
recipients.  The  intended  communicators  are  generally  unaware  that  this  is  occurring 
and  the  intent  of  the  unauthorized  listener  varies.  The  interceptor  may  remain  pas¬ 
sive  and  simply  “listen”  with  the  intent  of  monitoring,  recording,  analyzing,  etc., 
the  communication  activity.  This  type  of  passive  activity  is  very  difficult  to  detect. 
In  other  cases,  the  interceptor  may  become  active  and  “join”  in  the  communication 
activity.  This  may  take  the  form  of  “spoofing”  or  “man-in-the-middle”  type  attacks 
whereby  an  identity  compromise  occurs  and  the  unintended  party  is  able  to  freely 
inject  traffic  into  the  system.  This  activity  is  generally  detectable  given  that  inter- 
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ceptor  RF  emissions  are  present.  To  mitigate  this  activity,  there  is  a  pressing  need  to 
improve  both  pre-attack  security  and  post-attack  digital  forensics. 

1.1.2  Technical  Motivation. 

1.1. 2.1  PHY  Layer  Network  Security.  Much  research  has  focused  on 
traditional  bit-level  algorithmic  approaches  to  improve  network  security  and  mitigate 
spoofing.  More  recently  consideration  has  been  given  to  detecting  and  mitigating 
spoofing  near  or  at  the  bottom  of  the  Open  Systems  Interconnection  (OSI)  network 
stack.  One  such  work  includes  the  addition  of  a  “lightweight  security  layer”  hosted 
within  the  Medium  Access  Control  (MAC)  layer  to  detect  spoofing  and  anomalous 
traffic  [35].  Other  recent  efforts  have  focused  on  Physical  (PHY)  layer  implementa¬ 
tions  with  a  goal  of  exploiting  RF  characteristics  (radio  and  environmental)  that  are 
difficult  to  mimic,  thus  minimizing  the  opportunity  for  spoofing.  Two  such  efforts 
investigated  the  use  of  Received  Signal  Strength  (RSS)  as  a  means  for  detecting  the 
presence  of  a  spoofing  node  [5,53].  Although  related  in  their  use  of  RSS,  it  is  not 
entirely  clear  that  these  works  are  comparable  one-to-one  given  that  the  experiments 
were  conducted  using  different  hardware  being  operated  in  different  physical  environ¬ 
ments.  Under  dissimilar  conditions  such  as  these,  it  is  expected  that  statistics  of  the 
power-based  RSS  metric  would  vary.  This  variation  is  not  unique  to  wireless  com¬ 
munications  and  is  encountered  in  other  applications  employing  power-based  metrics, 
especially  when  all  system  and  environmental  interactions  are  accounted  for  (antenna 
patterns,  multipath,  background  noise,  etc.). 

The  authors  in  [53]  introduce  RF  fingerprinting  in  [24]  as  an  alternative  tech¬ 
nique  to  detect  and  mitigate  spoofing  through  PHY  layer  mechanisms.  However, 
they  readily  dismiss  this  alternative  for  “scale”  reasons.  Assuming  this  conclusion  is 
based  on  the  fact  that  RSS  is  currently  supported  and  provided  with  most  manufac¬ 
tured  devices,  the  authors’  position  is  supportable.  This  is  particularly  true  when 
constraining  PHY  layer  anti-spooffilg  mechanisms  to  reside  on  PC-sized  cards.  How¬ 
ever,  there  may  be  applications  where  the  size  constraints  are  much  more  relaxed 
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and  RF  fingerprinting  becomes  a  viable  alternative.  Consistent  with  related  work 
in  [54,55],  these  applications  were  addressed  through  the  fundamental  research  goal 
that  involved  demonstrating  radar- like  Specific  Emitter  Identification  (SEI)  capability 
similar  to  what  is  used  to  distinguish  between  radar  emitters  [6,9,11,12,34,40,46,60]. 

1.1. 2.2  Specific  Emitter  Identification.  Radar-based  SEI  research 
spans  nearly  twenty  years  and  has  considered  both  conventional  and  non-conventional 
parameters  for  identification.  Conventional  radar  parameters  generally  include  those 
which  are  based  on  intentional  modulation  which  may  be  applied  across  multiple 
pulses  (inter-pulse  modulation)  or  within  a  given  pulse  (intra-pulse  modulation). 
These  modulations  are  introduced  to  improve  some  aspect  of  overall  radar  perfor¬ 
mance  (tracking  accuracy,  ambiguity  resolution,  clutter  suppression,  etc.).  There  are 
other  unintentional  modulations  that  may  be  induced  by  the  hardware  used  to  im¬ 
plement  the  system  [30,34],  These  unintentional  modulations  may  result  from  any 
number  of  hardware  issues,  including  poor  system  design  (device  incompatibility), 
improper  operation  (over/under  voltage),  physical  device  limitations  (operating  tem¬ 
perature  range),  etc.  When  viewed  at  the  waveform  level,  many  of  these  features 
are  similar  to  what  currently  exist  in  modern  wireless  communication  systems  that 
typically  transmit  burst-like  waveforms  representing  various  forms  of  digital  informa¬ 
tion  (symbols,  bits,  packets,  etc.).  Communication  researchers  have  recognized  these 
similarities  and  have  begun  to  address  the  question:  “Can  existing  SEI  methods 
be  employed  with  wireless  communication  signals  to  achieve  radar-like  SEI 
capability?” 

1.1. 2.3  RF  Fingerprinting.  The  task  of  automatically  detecting,  iden¬ 
tifying  and  locating  commercial  RF  communication  devices  remains  a  challenging 
technical  problem.  The  work  presented  here  addresses  four  main  aspects  of  this 
problem,  including:  1)  the  selection  and  generation  of  fundamental  signal  charac¬ 
teristics  (amplitude,  phase,  and/or  frequency),  2)  the  feasibility  and  repeatability 
of  detecting  and  locating  the  start  of  a  burst  using  selected  waveform  feature(s) 
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amidst  channel  noise,  3)  the  identification  and  robust  extraction  of  distinguishable 
fingerprints-features  that  uniquely  characterize  the  unintentional  modulation  of  a 
device,  and  4)  the  performance  of  signal  classification  under  varying  channel  con¬ 
ditions  and  Signal-to-Noise  Ratio  (SNR).  The  ultimate  goal  is  to  demonstrate  an 
end-to-end  process  to  accurately  classify  commercially-available  RF  communication 
devices  using  signal  features  extracted  from  collected  emissions.  Relevant  research  in 
wireless  network  security  and  RF  fingerprinting  suggests  that  information  in  funda¬ 
mental  emissions  and  unintentionally  modulated  regions  provides  the  most  effective 
means  for  identifying  transmitters.  Collectively,  related  works  in  RF  fingerprint¬ 
ing,  electromagnetic  signatures,  intrapulse  modulation,  and  unintentional  modula¬ 
tion  [11,23,24,30,34,51,64,66,68],  form  a  solid  basis  for  developing  techniques  that 
may  be  applicable  to  commercial  communication  devices.  If  the  inherent  RF  finger¬ 
prints  are  repeatedly  extractable  and  sufficiently  unique,  they  are  potentially  useful 
for  determining  the  specific  make,  model,  and/or  serial  number  of  a  given  device. 

Previous  work  highlighted  signal  structure  uniqueness  and  attributed  inter¬ 
device  differences  to  various  manufacturing,  aging,  and  environmental  factors  [68]. 
While  several  processing  steps  are  required  to  effectively  exploit  the  unique  RF  fin¬ 
gerprints,  burst  location  is  arguably  the  most  important  [23,66].  In  this  context,  burst 
location  includes  determining  both  the  burst  start  time  and  the  subsequent  signal  re¬ 
gion^)  from  which  fingerprints  are  extracted.  Burst  detection,  burst  start  location 
and  signal  region(s)  selection  for  fingerprint  extraction  are  all  important  given  that 
improper  determination  of  the  burst  start  location  and  imprudent  selection  the  signal 
region(s)  can  adversely  bias  processing  in  favor  of  channel  noise  or  undesired  signal 
features  [68]. 

The  most  relevant  published  results  to  date  for  this  research  are  found  in  [54,55] 
and  is  based  on  experimentally  collected  802. 11A  signals.  As  with  this  previous 
work,  the  choice  of  using  Orthogonal  Frequency  Division  Multiplexing  (OFDM)-based 
signals  for  RF  fingerprinting  demonstration  was  driven  by  two  factors:  1)  consistency 
with  previous  related  802. 11A  work  that  has  been  extensively  published  [42,54,55, 
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63,67],  and  2)  the  continued  emergence  of  OFDM-based  signals  as  envisioned  for  4G 
Software  Defined  Radio  (SDR)  and  Cognitive  Radio  (CR)  communications  [21,26,48, 
72],  While  the  fundamental  fingerprinting  and  classification  techniques  in  this  work 
are  believed  to  be  broadly  applicable  to  other  signal  types,  the  challenges  posed  by 
OFDM-based  signals  is  of  near-term  interest. 

1.1. 2. 4  Signal  Denoising.  In  some  applications  the  desired  level  of 
performance  cannot  be  achieved  due  to  inherent  noise  contributions  in  the  environ¬ 
ment.  The  effective  mitigation  of  such  adverse  noise  effects  has  been  demonstrated 
in  numerous  applications  by  “denoising”  the  signal  of  interest  prior  to  processing 
to  remove  undesired  noise  contributions.  This  can  be  accomplished  using  a  Dis¬ 
crete  Wavelet  Transform  (DWT)  by  exploiting  differences  in  the  distribution  of  signal 
burst  energy  and  the  Additive  White  Gaussian  Noise  (AWGN)  in  which  it  is  embed¬ 
ded  [4,7,8,14,15,17-19,44,61].  The  common  approach  to  wavelet  denoising  includes: 
1)  transforming  the  input  signal  with  the  desired  transform, 2)  comparing  coefficient 
magnitudes  with  a  pre-defined  threshold,  3)  zeroing-out  all  coefficients  having  mag¬ 
nitudes  less  than  the  threshold  while  retaining  those  above  the  threshold,  4)  inverse 
transforming  the  thresholded  set  of  coefficients,  and  5)  processing  the  resultant  de- 
noised  signal. 

One  distinct  disadvantage  of  the  DWT  is  the  lack  of  shift  invariance,  i.e.,  for 
a  given  time  shift  in  the  input  signal  the  transformation  yields  a  different  set  of  co¬ 
efficients.  For  burst  detection,  this  problem  has  the  consequence  of  complicating  the 
computation  of  reasonable  thresholds  for  signal  denoising.  One  shift  invariant  (when 
properly  implemented)  alternative  is  the  Short  Time  Fourier  Transform  (STFT)  [13]. 
The  STFT  is  a  Fast  Fourier  Transform  (FFT)  done  over  a  series  of  short  contiguous 
time  intervals  spanning  the  signal  of  interest  [47].  Ideally,  the  intervals  are  short 
enough  to  maintain  piece-wise  stationarity  across  the  signal  while  at  the  same  time 
long  enough  to  capture  sufficient  spectral  energy.  This  trade-off  represents  a  compro¬ 
mise  between  achievable  time  resolution  (better  with  a  shorter  interval)  and  achievable 
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frequency  resolution  (better  with  a  longer  interval)-the  Heisenberg  inequality  [44], 
One  drawback  is  that  for  a  given  STFT  interval  length,  which  generally  remains  fixed 
throughout  the  signal  duration,  the  resolution  in  both  time  and  frequency  is  uniform 


across  the  domains.  This  is  illustrated  in  Figure  |T7I  )(a)  using  Heisenberg  uncertainty 
boxes. 


The  DWT  achieves  non-uniform  multi-resolution  capability  by  effectively  scal¬ 
ing  the  time  interval  inversely  proportional  to  frequency  such  that  a  relatively  narrow 
interval  is  used  to  capture  high  frequency  content.  This  is  illustrated  in  Figure [TTl|fb)| 
which  shows  representative  Heisenberg  uncertainty  boxes  for  an  arbitrary  DWT.  As 
indicated,  the  higher  frequency  content  regions  have  higher  time  resolution  (narrower 
box  widths  across  time)  and  lower  frequency  content  regions  have  lower  time  reso¬ 
lution  (wider  box  widths  across  time).  This  type  of  multi-resolution  time- frequency 
decomposition  works  best  if  the  signal  is  composed  of  high  frequency  components  of 
short  duration  plus  low  frequency  components  of  long  duration,  characteristics  which 
most  signals  possess  [47]. 

An  alternative  wavelet  transform  that  possesses  both  the  DWT’s  multi-resolution 
capability  and  the  STFT's  shift  invariance  is  the  Dual-Tree  Complex  Wavelet  Trans¬ 
form  (DT-CWT)  [2,50].  The  DT-CWT  is  a  DWT  extension  that  is  “nearly  shift- 
invariant,”  i.e. ,  the  DT-CWT  coefficients  are  independent  of  time  domain  shift  and 
more  strongly  dependent  on  inter-scale  and  intra-scale  neighborhoods  [50].  Further¬ 
more,  the  DT-CWT  magnitude  response  exhibits  reduced  ringing  in  the  wavelet  do¬ 
main  due  to  high-frequency  noise  and  sharp  discontinuities,  which  makes  the  denoising 
process  more  reliable  by  ensuring  consistent  threshold  calculations  [50]. 


1.2  Research  Contributions 

Table  11.11  provides  a  list  of  various  Technical  Areas  (concepts,  techniques,  at¬ 
tributes,  metrics,  etc.)  and  the  relational  mapping  between  Previous  related  work 
and  the  Current  research  presented  in  this  dissertation.  As  summarized  in  the  fol- 
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Frequency 


(a)  Short  Time  Fourier  Transform  (STFT) 


(b)  Discrete  Wavelet  Transform  (DWT) 


Figure  1.1:  Tiling  of  Heisenberg  uncertainty  boxes  in  the  time-frequency  plane  for 
STFT  and  DWT  decompositions.  The  width  and  height  of  a  given  box  is  related  to 
its  time  and  frequency  resolution,  respectively  [47]. 


lowing  subsections,  there  have  been  contributions  made  to  each  of  the  technical  areas 
identified  in  the  first  five  rows  of  the  table. 

1.2.1  Performance  Criteria.  Experimental  setup  and  execution  can  differ 
from  one  research  activity  to  another,  even  when  using  identical  or  similar  equipment 
and  processes.  This  can  make  direct  comparison  of  new  results  with  previous  results 
difficult  and  care  must  be  taken  to  ensure  that  1)  previous  contributions  are  fairly 
represented  and  2)  new  contributions  are  sufficiently  supported-this  is  the  case  for 
work  presented  here.  For  all  previous  AFIT-based  works  referenced  in  Table  11.11 
there  were  no  firm  performance  goals  or  criteria  in  place  at  the  time  the  work  was 
conducted.  Rather,  proof-of-concept  demonstration  was  the  main  objective  and  “As 
Achieved”  performance  was  reported  as  noted  in  Table  11.11 

While  as  achieved  results  are  presented  in  this  document  as  well,  and  based  on 
many  combinations  of  parameters  and  parameter  values,  specific  performance  criteria 
was  introduced  to  help  highlight  performance  differences  (poorer  and  better)  across 
the  numerous  scenarios  considered.  The  “Reasonable”  criteria  used  here  and  shown 
in  Table  fTTl is  somewhat  arbitrary  and  based  on  achieving  80%  or  better  classification 
accuracy  at  SNR  <  20  dB.  Using  this  reasonable  operating  point  of  80%  classification 
accuracy,  performance  comparisons  are  made  throughout  Chapter  4  based  on  the 
“gain”  provided  by  Wavelet  Domain  (WD)  techniques  relative  to  what  is  provided  by 
Time  Domain  (TD)  techniques.  This  gain  is  defined  here  as  the  reduction  in  required 
SNR,  in  dB,  for  the  WD  fingerprinting  technique  to  achieve  the  same  classification 
performance  as  the  TD  fingerprinting  technique. 

1.2.2  TD  Fingerprint  Classification.  Prior  to  assessing  WD  classification 
performance,  it  was  necessary  to  replicate  TD  results  in  [54]  to  form  a  baseline  for 
comparison.  Upon  replicating  these  earlier  results,  it  became  evident  that  the  post- 
collection  filter  bandwidth  BWpc  was  a  very  important  parameter  and  that  all  earlier 
TD  results  were  based  on  using  a  fixed  value.  While  the  fixed  bandwidth  approach  was 
sound  and  the  selected  bandwidth  was  reasonably  based  on  sound  engineering  prac- 


Table  1.1:  Relational  mapping  between  Technical  Areas  in  Previous  related  work 
and  Current  research  contributions. 


Technical  Area  Previous  Current 


Addressed 

Ref# 

Addressed 

Ref# 

TD  Fingerprinting 

X 

[23,24,54,55,68] 

X 

[31-33] 

WD  Fingerprinting 

X 

[32] 

SNR  Sensitivity 

X 

[54, 55] 

X 

[31-33] 

Burst  Detection 

X 

[31,33] 

Signal  Type  /  Modulation 

802.11A  /  OFDM 

X 

[20,42,54,55,63,67] 

X 

[31-33] 

802.11B  /  DSSS 

X 

[24,66,68] 

802.11G  /  OFDM 

X 

GSM  /  GMSK 

X 

[3] 

Bluetooth  /  GFSK 

X 

[23,71] 

Instantaneous  Signal  Characteristics 

Amplitude 

X 

[20,24,42,51] 

[63,64,66-68] 

X 

[31-33] 

Phase 

X 

[20,23,24] 

X 

[31-33] 

Frequency 

X 

[20,24] 

X 

[31—33] 

RF  Fingerprint  Features  and  Metrics 

Std  Deviation 

X 

[20,24] 

Variance 

X 

[20,23] 

X 

[31-33] 

Skewness 

X 

[54] 

X 

[31—33] 

Kurtosis 

X 

[54] 

X 

[31-33] 

Classification  Method  and  Performance  Criteria 


Bayesian  MDA/ML 

X 

[20,54] 

X 

[31-33] 

“As  Achieved” 

X 

[23,24,51,54,68] 

“Reasonable” 

X 

[32] 
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tices,  the  earlier  works  provided  no  bandwidth  sensitivity  analysis.  This  analysis  was 
subsequently  carried  out  under  this  research  and  a  bandwidth  of  BWpc  =  7.7  MHz 
was  used  for  generating  all  comparative  TD  and  WD  results.  This  particular  value 
enabled  comparison  of  both  techniques  at  their  best  overall  performance  levels,  with 
TD  having  an  approximate  2%  advantage  in  device  classification  at  higher  SNRs- 
an  advantage  that  rapidly  diminishes  at  lower  SNRs  that  are  more  consistent  with 
operational  environments  [31-33]. 

1.2.3  WD  Fingerprint  Classification.  Relative  to  TD  fingerprint  classifica¬ 
tion,  enhanced  fingerprint  classification  is  demonstrated  here  using  improved  finger¬ 
print  features.  Specifically,  this  work  represents  the  first  application  of  a  DT-CWT 
decomposition  to  enhance  features  of  statistical  RF  fingerprints.  Considerable  per¬ 
formance  improvement  or  gain  is  realized  using  the  enhanced  WD  feature  set  with 
identical  post-collection  filter  bandwidth  and  Multiple  Discriminant  Analysis/ Maxi¬ 
mum  Likelihood  (MDA/ML)  processing  [32], 

1.2.4  SNR  Sensitivity  Analysis.  With  the  exception  of  results  generated 
under  this  research  and  documented  in  [31-33],  a  majority  of  the  works  cited  in 
Section  11.1.21  lack  any  form  of  sensitivity  analysis  in  terms  of  assessing  burst  detec¬ 
tion  and/or  fingerprint  classification  performance  under  varying  channel  noise  con¬ 
ditions.  The  two  exceptions  are  the  most  recent  related  works  in  [54,55].  Noise 
sensitivity  analysis  is  imperative  for  determining  acceptable  SNR  levels  for  achiev¬ 
ing  consistent  and  reliable  classification  results.  Classification  sensitivity  to  channel 
noise  and  burst-to-burst  detection  variability  has  been  analyzed  using  experimen¬ 
tally  collected  802. 11 A  signals  in  [33].  With  respect  to  burst  location  estimation, 
both  Fractal-Bayesian  Step  Change  Detector  (BSCD)  and  Traditional  Variance  Tra¬ 
jectory  (VT)  techniques  provided  results  that  were  consistent  with  “perfect”  burst 
location  (a  start  location  based  on  visual  inspection  of  each  collected  burst)  at  higher 
SNRs  (10  <  SNR  <  30  dB).  However,  performance  for  both  techniques  diverged  at 
lower  SNRs  (—3  <  SNR  <  10  dB)  [33].  With  respect  to  the  burst  location  esti- 
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mation  error  impact  to  classification  performance,  the  Traditional  VT  technique  was 
consistent  with  perfect  estimation  for  6  <  SNR  <  30  dB  but  underperformed  for 
—3  <  SNR  <  6  dB.  Traditional  VT  also  provided  considerable  improvement  when 
compared  with  the  Fractal-BSCD  technique  at  lower  SNRs  (—3  <  SNR  <  18  dB), 
i.e. ,  for  a  given  classification  accuracy  in  the  range  of  50%-80%  the  required  SNR  for 
Traditional  VT  is  3-6  dB  lower  than  what  is  required  for  Fractal-BSCD.  This  short¬ 
fall  provided  an  impetus  for  subsequent  burst  detection  research  aimed  at  improved 
performance  at  low  SNRs  [31]. 

1.2.5  Burst  Detection  at  Lower  SNR.  As  published  in  [31]  and  presented  in 
this  dissertation,  signal  denoising  with  the  DT-CWT  prior  to  Traditional  VT  burst 
detection  (introduced  here  as  Denoised  VT  processing)  is  more  effective  and  provides 
more  robust  burst  detection  and  location  at  lower  SNRs  (—3  <  SNR  <  10  dB). 
Relative  to  results  for  perfect  burst  detection  and  location,  the  Denoised  VT  pro¬ 
cess  achieves  nearly  34%  of  the  available  performance  improvement-when  used  with 
MDA/ML  processing,  there  is  little  more  to  be  gained  in  overall  classification  perfor¬ 
mance  by  improving  burst  detection  and  location  accuracy  [31]. 

1.3  Dissertation  Overview 

This  document  is  divided  into  five  chapters  and  contains  one  appendix.  Chap¬ 
ter  2  presents  relevant  technical  background  information  on  major  concepts  and  tech¬ 
niques  used  to  conduct  the  research.  Sufficient  technical  detail  is  presented  such  that 
the  fundamental  research  approach  is  repeatable  and  the  key  contributions  are  verifi¬ 
able.  The  major  concepts  and  techniques  are  presented  as  functionally  implemented 
in  the  overall  demonstration  process. 

Chapter  3  provides  the  overall  demonstration  process  used  for  generating  results 
and  conducting  analysis.  A  detailed  description  is  included  for  both  the  “Signal  Col¬ 
lection”  hardware  and  “Post-Collection  Processing”  software  processes.  The  primary 
hardware  used  for  signal  collection  was  AFIT’s  RF  Signal  Intercept  and  Collection 
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System  (RFSICS)  with  subsequent  data  processing  accomplished  exclusively  in  a 
MATLAB®  environment. 

Chapter  4  provides  modeling,  simulation  and  analysis  results  that  were  gener¬ 
ated  using  the  processes  detailed  in  Chapter  3.  The  research  involved  hundreds  of 
simulations,  each  requiring  tens  of  hours  of  processing  time  in  some  cases.  For  brevity 
and  to  ensure  succinctness,  only  a  subset  of  representative  results  are  presented  from 
selected  scenarios  to  fully  support  key  research  findings  and  contributions. 

Chapter  5  concludes  the  main  document  by  providing  an  overall  summary  of 
research  activities,  a  summary  of  key  findings,  and  recommendations  for  subsequent 
research.  This  is  followed  by  an  appendix  that  provides  some  of  the  developmental 
MATLAB®  code  used  to  support  the  research. 
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II.  Background 

This  chapter  presents  relevant  technical  background  information  on  major  concepts 
and  techniques  used  to  conduct  the  research.  The  material  here  supports  subsequent 
material  presented  in  the  methodology,  results  and  conclusion  chapters  of  the  docu¬ 
ment.  This  chapter  is  not  presented  as  a  complete  tutorial,  but  rather,  intended  to 
provide  sufficient  detail  such  that  the  fundamental  research  approach  is  repeatable 
and  the  key  contributions  are  verifiable.  For  convenience,  the  major  concepts  and 
techniques  are  presented  as  functionally  implemented  in  the  overall  demonstration 
process.  Burst  Detection  and  Location  is  first  presented  in  Section  l2Tl  which  provides 
details  on  the  two  specific  techniques  considered,  including  the  Fractal-Bayesian  Step 
Change  Detector  (Fractal-BSCD)  in  Section  f2. 1.1 1  and  the  Traditional  Variance  Tra¬ 
jectory  (Traditional  VT)  technique  in  Section  [2.1.21  Lastly,  the  Dual- Tree  Complex 
Wavelet  Transform  (DT-CWT)  is  presented  in  Section  1231 

2.1  Burst  Detection  and  Location 

Discriminating  a  burst-like  signal  response  from  background  noise  can  be  a 
difficult  task  as  many  burst  responses  can  appear  noise-like.  In  some  respects,  de¬ 
tecting  a  burst  response  is  akin  to  separating  noise  from  noise  [52],  Related  research 
has  focused  on  exploiting  two  different  properties  to  discriminate  between  signal  and 
background  noise  contributions,  including  inherent  signal  structure  and  instantaneous 
signal  characteristics.  Inherent  signal  structure  has  been  successfully  exploited  using 
a  Fractal-Bayesian  Step  Change  Detector  (Fractal-BSCD)  while  instantaneous  signal 
characteristics  have  been  exploited  using  Traditional  Variance  Trajectory  (VT).  The 
details  for  these  approaches  are  provided  in  the  following  subsections. 

2.1.1  Fractal-Bayesian  Step  Change  Detector.  The  Fractal-Bayesian  Step 
Change  Detector  (Fractal-BSCD)  has  been  used  to  exploit  inherent  signal  structure  to 
discriminate  between  signal  and  channel  background  responses.  As  time  progresses, 
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random  Additive  White  Gaussian  Noise  (AWGN)  exhibits  no  structure  amongst  sam¬ 
ples  while  a  deterministic  signal  does.  When  a  time  series  transitions  from  a  region 
containing  only  noise  to  a  region  containing  both  noise  and  signal  its  inherent  struc¬ 
ture  changes.  This  change  can  be  detected  using  fractal  dimensions.  However,  a  burst 
response  in  such  a  region  is  non-stationary  and  therefore,  not  a  pure  fractal,  i.e. ,  its 
fractality  is  a  function  of  time  and  thus  it  cannot  be  self-similar  [52],  Otherwise,  the 
calculated  fractal  dimension  would  yield  the  same  value  regardless  of  the  signal  time 
and  duration  used,  which  does  not  describe  a  non-stationary  signal.  Yet  on  a  smaller 
scale,  a  transient  can  have  local  stationary  fractality  and  can  be  modeled  as  a  series 
of  piece-wise  fractals  through  multi-fractality  analysis.  The  local  fractal  dimensions 
are  calculated  using  a  sliding  window  [52], 

It  has  been  demonstrated  that  burst  start  location  can  be  accomplished  using 
the  fractal  dimension  [64]  measure  followed  by  a  Bayesian  Step  Change  Detector  [42, 
63,64,67].  This  process  is  denoted  here  as  Fractal-BSCD.  The  fractal  derivation  can 
be  found  in  [27]  and  can  be  calculated  using  the  following  Higuchi  method.  Given  a 
windowed  data  time  series  (A"(l),  X(2),  ...,  X(Nx)},  the  curve  length  is  defined  as: 


(2.1) 


nl 

X  =  \x  {m  +  ik)  -  X  (m  +  [i  -  1)  k)  |  , 

i= 1 

where  Nl  =  \_(NX  —  m)/k\ ,  [_*J  is  the  floor  operator,  k  is  the  interval  index  number, 
and  me  [1,  k\  is  the  start  time  index  number. 

The  average  of  Lm  (k)  over  m  is  denoted  as  (L  (k))  and  defines  the  curve  length 
for  time  interval  k.  By  varying  k  over  [1,  kmax]  and  plotting  ( L(k ))  versus  k  on  a 
log- log  scale,  the  data  ideally  forms  a  straight  line,  with  a  proper  selection  of  kmax . 
The  fractal  dimension  d  is  defined  as  the  negative  of  the  line  slope,  which  can  be 
calculated  using  a  least  squares  method.  Furthermore,  kmax  is  empirically  chosen.  If 
it  is  too  large,  the  data  plotted  on  the  log-log  scale  will  not  be  linear.  If  it  is  too 
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small,  there  will  not  be  enough  data  points  for  an  accurate  linear  fit.  For  this  work, 
a  value  of  /cmax  =  10  is  chosen  for  all  fractal  calculations. 

Using  the  fractal  dimension  vector  d  formed  across  all  data  windows,  BSCD 
is  applied  to  determine  the  a-posteriori  probability  that  a  given  fractal  dimension 
dm  E  d  represents  the  data  change  point  corresponding  to  the  burst  start.  The 
a-posteriori  Probability  Distribution  Function  (PDF)  for  m  given  d  is  [41] 
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where  Np  is  the  length  of  d,  [_•]  is  the  floor  operator,  /  denotes  prior  information, 
and  m  is  the  potential  change  point  being  evaluated.  The  value  of  m  corresponding 
to  max[p({m\  |d,  /)]  establishes  the  burst  start  sample  number.  Representative  re¬ 
sponses  for  Fractal-BSCD  processing  are  shown  in  Figure  EH]  where  the  circled  region 
highlights  the  burst  start  location  at  t  —  0.  As  illustrated  in  the  bottom  a-posteriori 
PDF  response  there  is  a  distinct  peak  that  corresponds  to  the  burst  start  time. 


Work  in  [23,  68]  shows  that  abrupt,  non-gradual  feature  changes  are  impor¬ 
tant  for  the  Fractal-BSCD  process  to  work  effectively.  Signals  having  more  gradual 
ramp-like  versus  impulse-like  responses  are  problematic  and  require  alternate  meth¬ 
ods  of  detection.  Similar  BSCD-based  methods  have  been  considered  to  address  the 
increased  challenge,  e.g.,  Bayesian  Ramp  Change  Detection  [63,67].  However,  as  de¬ 
tailed  in  the  next  section  there  are  alternatives  to  BSCD-based  methods  that  have 
proven  effective  as  well. 


2.1.2  Traditional  Variance  Trajectory.  The  Traditional  Variance  Trajectory 
(Traditional  VT)  alternative  to  burst  detection  exploits  instantaneous  signal  charac¬ 
teristics  to  discriminate  between  signal  and  channel  background  responses.  While 
the  Traditional  VT  process  can  be  applied  to  any  arbitrary  sequence  of  data,  it  has 
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Figure  2.1:  Representative  responses  for  Fractal-BSCD  processing:  (Top)  Instan¬ 
taneous  Amplitude,  (Middle)  Fractal  d,  and  (Bottom)  A-Posteriori  PDF.  The  circled 
region  highlights  the  burst  start  location  at  t  —  0. 

previously  been  used  for  burst  detection  with  both  instantaneous  phase  [23]  and  in¬ 
stantaneous  amplitude  characteristics  [55].  Given  an  arbitrary  input  sequence,  the 
Traditional  VT  process  consists  of  1)  dividing  the  input  sequence  into  sequential  sub¬ 
sequences,  or  windows  of  data,  which  may  or  may  not  overlap,  2)  calculating  the 
variance  over  each  window  of  data,  and  3)  forming  the  “trajectory”  sequence  as  the 
difference  between  consecutive  window  variances.  Given  arbitrary  sequence  {x(k)}, 
k  =  1,2,...,  Nx,  the  variance  trajectory  of  {x{k)}  is  denoted  as  the  sequence  {VTx(i)} 
where  the  ith  element  is  given  by  [55] 


VTx(i)  =  \Wx(i)-Wx{i  +  1)|  , 


i  1,2,...,  Lw  1  , 


(2.3) 
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l  +  (m - l^Ns^Nw 

wx{m)  =  —  ix(k)  -  >  (2-4) 

“  fc=l+(m-l)JVs 

77i  1,2,...,  Lw  , 

where  AC,  is  the  window  extent,  and  Ns  is  the  number  of  samples  the  window  advances 
between  sequential  calculations.  The  [iw  factor  in  (12.41  is  the  sample  mean  of  { xw(k )} 
which  is  the  subsequence  of  consecutive  elements  from  {x(k)}  contained  in  window 

w. 

Figure  [2721  shows  representative  responses  for  Traditional  VT  processing  where 
the  top  plot  is  the  magnitude  response  of  {x{k)}  and  the  other  two  plots  are  the 
corresponding  responses  for  Traditional  VT  at  SNR  =  40  dB  and  SNR  =  0  dB.  The 
circled  region  highlights  the  burst  start  location  at  t  =  0.  As  seen  in  the  SNR  = 
40  dB  response,  there  is  a  distinct  peak  corresponding  to  the  burst  start  time  near 
t  =  0.  The  sensitivity  of  Traditional  VT  processing  to  SNR  variation  is  evident  in  the 
SNR  =  0  dB  response  where  the  peak  response  near  the  burst  start  time  is  virtually 
indistinguishable  from  earlier  ( t  <  0)  peaks.  As  used  here  and  in  other  previous 
work  with  instantaneous  signal  characteristics,  the  degradation  of  Traditional  VT 
performance  at  lower  SNRs  directly  impacts  burst  detection  and  location  error  and 
subsequent  classification  performance. 

2.2  RF  Fingerprint  Classification 

There  has  been  considerable  work  in  previous  years  involving  the  exploitation 
of  RF  signal  characteristics  to  classify  signals  and  identify  the  devices  producing 
them  [23,54,55,64,66,68].  Collectively,  these  works  embody  the  held  of  RF  Fingerprint 
Classification  which  fundamentally  requires  two  processes,  including:  1)  fingerprint 
generation  and  2)  fingerprint  classification.  Fingerprint  generation  requires  the  se¬ 
lection  and  extraction  of  features  that  enable  signal/device  discrimination.  Desirable 
properties  of  the  selected  feature  set  include:  1)  reduced  dimensionality  to  minimize 
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Time  (|isec) 


Figure  2.2:  Representative  responses  for  VT  processing:  (Top)  Instantaneous  Am¬ 
plitude,  (Middle)  SNR  =  40  dB,  and  (Bottom)  SNR  =  0  dB.  The  circled  region 
highlights  the  burst  start  location  at  t  —  0. 


processing  and  storage  requirements,  2)  intra-device  repeatability,  and  3)  inter-device 
uniqueness.  For  this  work,  the  classification  features  are  statistics  of  instantaneous 
signal  characteristics  per  the  details  provided  in  Section  12.2.11  and  Section  12.2.21  The 
resultant  RF  Statistical  Fingerprints  are  then  used  for  signal/device  classification  per 
the  details  provided  in  Section  12.2.31 


2.2.1  Instantaneous  Signal  Characteristics.  While  there  are  many  signal 
characteristics  that  could  be  used  for  device  identification  (instantaneous  responses, 
peak  responses,  average  responses,  amplitude,  phase,  frequency,  power,  etc.),  a  ma¬ 
jority  of  earlier  related  works  have  predominantly  focused  on  instantaneous  amplitude 
and  instantaneous  phase  characteristics  [23,64,66,68].  The  most  recent  research  has 
exploited  instantaneous  frequency  characteristics  as  well  [54, 55] .  As  adopted  for  con¬ 
sistency  with  these  previous  work,  the  following  development  of  instantaneous  signal 
characteristics  is  provided  for  completeness. 
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Samples  of  a  complex  time  domain  (TD)  signal  having  in-phase  and  quadrature 
components  of  ItdQi)  and  Qtd  ('»-),  respectively,  can  be  expressed  as 


std (n)  —  Itd{p)  +  jQtd(^)  ,  (2.5) 

and  have  corresponding  instantaneous  amplitude,  a(n),  instantaneous  phase,  4>(n), 
and  instantaneous  frequency,  f(n),  responses  are  given  by 

a  (n)  =  \J Itd  (n)  "h  Qtd  ( n )  >  (2-6) 


(f)  (n)  =  tan  1 


Qtd  (n) 
Itd  (n)  _ 


(2.7) 


1  4>(n)  —  (f>(n  —  1) 

27 r  An 


(2.8) 


In  practice,  each  characteristic  response  is  “centered”  (mean  removed)  to  re¬ 
move  collection  system  biases  that  may  unduly  influence  subsequent  processing.  The 
instantaneous  amplitude  and  frequency  responses  are  simply  centered  using 


ac{n)  =  a(n)  -  fxa  , 


(2.9) 


fc{n)  =  f(n )  -  a if  ,  (2.10) 

where  n  —  1,2,3,...,  Nm,  Nm  is  the  total  number  of  samples  in  the  sampled  signal, 
and  /ia  and  // /  are  amplitude  and  frequency  means  calculated  across  Nm  samples  of 
(12 .Oft  and  (12.81).  respectively. 

The  phase  centering  process  is  somewhat  more  involved  and  includes  removal  of 
a  linear  phase  component  prior  to  centering.  This  component  may  be  due  to  collection 
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receiver  coloration  or  result  from  inexact  frequency  estimation  during  post-collection 
down-conversion.  Given  the  phase  response  in  (12. 71).  the  non-linear  phase  response  is 
given  by 


4>ni(n)  =  0(n)  -  2nfj,f(n)At  , 


(2.11) 


where  /j,f  is  the  frequency  mean  used  in  (12. 10ft  and  At  is  the  time  sample  spacing.  As 
a  final  step,  the  mean  of  cpni  is  removed  to  yield  the  desired  centered  non-linear  phase 
which  is  given  by 


(ficnl (^)  (pnlijl)  > 


(2.12) 


where  /i^  is  the  mean  of  (j)ni{n)  in  (12. lip.  The  centering  of  signal  characteristics  in 
(I2.9p~02.12ft  is  consistent  with  previous  fingerprint  classification  work  that  successfully 
employed  similar  procedures  [54,55]. 

2.2.2  Statistical  Feature  Metrics.  Direct  use  of  signal  characteristics  such 


as  those  presented  in  Section  12.2.11  for  classification  features  can  be  prohibitive  in 


terms  of  data  storage  memory  requirements  and  computational  processing  time.  The 


computational  burden  can  be  eased  by  reducing  the  feature  dimensionality  used  for 


fingerprint  classification.  This  was  successfully  accomplished  in  previous  work  using 
inherent  statistical  behavior  of  the  signal  characteristics  vice  the  signal  characteristics 
themselves  [54,55].  As  adopted  from  this  earlier  work,  the  statistics  of  interest  here 
included  the  variance  (a2),  skewness  (7),  and  kurtosis  {k).  For  arbitrary  sequence 
{x(A;)},  k  —  1,2,  ...,NX,  these  statistics  are  defined  as  [36]: 


(2.13) 
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N, 


(2.14) 


Nx 

E  ix(k)  -  xf 

k= 1 _ 

1  f  N*  ,)V2  ’ 

JV* 

i  E  [*(*)  -  *]4 

^  ,  (2-15) 

{is 

where  x  is  the  sample  mean  of  (x(fc)}.  The  final  RF  statistical  fingerprints  are 
formed  by  calculating  these  statistics  for  the  appropriate  centered  instantaneous  signal 
characteristic(s)  in  Section  12.2.11  i.e.,  setting  {x(k)}  equal  to  (ac(n)}  with  elements 
from  (I2.9F  setting  {x(k)}  equal  to  {/c(u)}  with  elements  from  (12.1  OF  and/or  setting 
{x(fc)}  equal  to  {ficnifn)}  with  elements  from  (I2.12F 

2.2.3  MDA/ML  Classification.  While  many  different  techniques  have  been 
researched  and  are  available  for  classification,  they  all  employ  two  fundamental  pro¬ 
cesses:  training  and  classification.  That  is,  they  train  the  classifier  using  a  subset 
of  the  input  data  and  then  classify  using  the  remaining  data.  For  the  most  part, 
these  techniques  are  oblivious  to  what  the  input  data  actually  represents  and  their 
performance  is  predominantly  driven  by  the  statistical  behavior  of  the  data.  With 
regard  to  RF  fingerprint  classification,  there  has  been  little  novelty  in  developing  spe¬ 
cialized  classification  techniques  and  most  researchers  have  opted  for  well-established 
techniques.  The  predominant  techniques  of  choice  have  been  based  on  neural  net¬ 
works  [45,51,52,57-59,62,63,65],  with  some  limited  additional  work  based  on  Kalman 
filtering  and/or  a  Hotelling  statistic  [22,28]. 

Multiple  Discriminant  Analysis/Maximum  Likelihood  (MDA/ML)  classification 
has  emerged  as  a  viable  alternative  and  successfully  used  for  RF  fingerprint  classifi¬ 
cation  [54].  Multiple  Discriminant  Analysis  (MDA)  is  an  extension  of  Fisher’s  Linear 
Discriminant  (FLD)  process  for  more  than  two  classes  [16].  For  a  3-class  problem,  the 
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MDA  process  projects  higher-dimensional  data  onto  a  2-dimensional  “Fisher  plane” 
that  maximizes  inter-class  distances  while  simultaneously  minimizing  intra-class  dis¬ 
tances.  In  principle,  this  method  cannot  improve  classification  potential.  However,  it 
provides  good  class  separation  and  visualization  of  data  having  input  dimensionality 
greater  than  three.  Using  this  lower-dimensional  data,  decision  boundaries  calcu¬ 
lated  from  ML  distributions  are  determined  assuming  normally  distributed  input 
data,  equal  costs  and  uniform  prior  probabilities.  In  general,  to  discriminate  c  classes 
using  d- dimensional  input  data,  the  input  vector  x  is  linearly  projected  onto  a  (d  —  1)- 
dimensional  space  using 

y  =  WTx  ,  (2-16) 

where  y  is  the  vector  of  projected  values  and  W  is  a  d  x  (c  —  1)  projection  ma¬ 
trix.  Classification  is  performed  using  unknown  data  and  the  trained  2-dimensional 
decision  boundaries  calculated  from  ML  distributions.  The  process  classifies  each 
“unknown”  input  data  set  by  projecting  it  onto  the  trained  Fisher  plane  according  to 
(12.161).  Projected  points  falling  within  the  correct  region  are  correctly  classified  while 
those  falling  outside  the  correct  region  are  misclassified.  The  percentage  of  correct 
classification  is  determined  based  on  the  total  number  of  unknown  trials.  A  more 
complete  description  of  the  MDA/ML  process  is  provided  in  [10]. 

2.3  Dual-Tree  Complex  Wavelet  Transform 

Device  classification  can  be  performed  using  a  Discrete  Wavelet  Transform 
(DWT),  with  one  popular  method  using  a  subset  of  the  largest  DWT  coefficient  mag¬ 
nitudes  as  the  classification  features  [44],  As  mentioned  in  Sectionfl. 1.2.41  one  distinct 
disadvantage  of  DWT-based  approaches  is  that  the  DWT  is  not  shift  invariant.  As 
with  signal  denoising,  this  presents  a  problem  for  RF  fingerprinting  applications  given 
that  robust  classification  performance  relies  on  the  fingerprint  features  being  unique, 
repeatable  and  stable.  These  properties  cannot  be  assured  if  the  underlying  features 
(DWT  coefficients)  vary  dramatically  throughout  the  processing  interval  of  interest. 
For  example,  variation  in  burst  detection  and  start  location  error  generally  translates 
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Figure  2.3:  Four  Stage  Dual- Tree  Complex  Wavelet  Transform  (DT-CWT)  [2]. 

to  greater  variation  in  fingerprint  features.  To  address  the  lack  of  shift  invariance  in 
DWT  processing,  a  Dual- Tree  Complex  Wavelet  Transform  (DT-CWT)  is  considered. 

The  DT-CWT  is  a  DWT  extension  that  is  “nearly  shift-invariant,”  i.e. ,  the 
DT-CWT  coefficients  are  independent  of  time  domain  shift  and  more  strongly  de¬ 
pendent  on  interscale  and  intrascale  neighborhoods  [50].  This  shift  invariance  has 
been  previously  exploited  to  improve  classification  performance  for  hyperspectral  im¬ 
ages  [38].  Furthermore,  the  DT-CWT  magnitude  response  exhibits  reduced  ringing 
that  is  generally  induced  by  high-frequency  noise  and  sharp  discontinuities  [50]. 

The  DT-CWT  is  commonly  implemented  using  two  real-valued  filter  banks. 
These  are  denoted  as  Treel  and  Tree 2  in  Figure  12.31  which  shows  one  common  ar¬ 
chitecture  for  DT-CWT  implementation  [2],  The  scaling  and  wavelet  functions  for 
Treel  are  symmetric  (even  functions)  while  Tree 2  has  scaling  and  wavelet  functions 
that  are  anti-symmetric  (odd  functions).  The  wavelet  and  scaling  functions,  ip(t)  and 
c t>(t )  respectively,  for  the  Treel  filter  bank  are  given  by  [2,50], 


W)  =  hi(n)(j)(2t  —  n)  , 

n 


(2.17) 
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(2.18) 


(j){t)  =  \Fl  h0(n)(j)(2t  -  n)  , 

n 

where  the  filter  coefficients  h\{n)  and  ho(n)  arc  implemented  directly  as  the  Analysis 
Filters  (AF)  given  in  [49]  (see  Section  [A.  511.  Ideally,  the  corresponding  functions  for 
the  Tree2  filter  bank  are  the  Hilbert  transforms  of  (12.171)  and  (12.181).  expressed  as 


-0  (t)  =  V2  22  ^i('«-)0/(2t  —  n)  , 

n 

(2.19) 

</>'(t)  =  y/2  *22  h'0(n)(f>,(2t  -  n )  , 

(2.20) 

n 


where  the  filter  coefficients  h\  (n)  and  h'0(n)  arc  implemented  directly  as  the  Analysis 
Filters  (AF)  given  in  [49]  (see  Section  fA.5l). 

As  shown  in  Figure  [2731  the  first  stage  filters  for  both  Treel  and  Tree2  have  dif¬ 
ferent  coefficients  when  compared  to  the  later  stage  filters  and  are  denoted  as  h^\n), 
h^\n),  h'2\n),  and  liQl\n)i  respectively.  The  first  stage  hlter  coefficients  are  im¬ 
plemented  directly  as  the  First  Analysis  Filters  (FAF)  given  in  [49]  (see  Section  fA. 51). 

For  real- valued  input  signals,  the  Treel  and  Tree2  hlter  banks  yield  real- valued 
wavelet  domain  (WD)  coefficients  representing  real  (. IlWD )  and  imaginary  (QlWD)  com¬ 
ponents  of  complex  coefficients  [50] .  These  components  can  be  functionally  combined 
in  a  form  similar  to  (12.51)  and  expressed  as 

sWD^n)  =  I\VD(n)  +  jQwoi.71)  ■  (2-21) 

Using  slWD(n )  elements  from  (12.21)).  the  sequence  {swr>(n)}  of  all  elements 
can  be  interpreted  as  what  may  be  called  a  “complex  sampled  WD  signal.”  Given 
the  similar  structure  of  this  WD  signal  and  the  TD  signal  in  (12.51),  WD  fingerprint 
classification  can  be  performed  using  the  process  in  Section  [2721  In  this  case,  the  WD 
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signal  in  (12.21)1  can  be  used  in  (12. 61)— (12 . 1 2D  to  generate  WD  signal  characteristics  and 
statistics  calculated  per  (12. 13)1— (12. 15)1  to  form  statistical  WD  fingerprints. 

2-4  Denoising 

Wavelet  transforms,  and  in  particular  the  DWT,  have  been  used  to  denoise 
signals  by  exploiting  differences  in  the  distribution  of  signal  and  embedded  noise 
contributions  in  the  wavelet  domain  [4,7,8,14,15,17-19,44,61].  In  the  case  of  an 
AWGN  channel,  the  noise  contribution  remains  Gaussian  in  the  wavelet  domain  and 
thus  uniformly  distributed  with  respect  to  scale  [43].  However,  burst  signal  contri¬ 
butions  are  non-uniformly  distributed  in  the  wavelet  domain  and  significant  signal 
content  is  generally  manifested  in  large  wavelet  coefficient  magnitudes.  Thus,  one 
common  approach  for  denoising  using  wavelets  involves  1)  transforming  the  time 
domain  signal  into  the  wavelet  domain,  2)  thresholding  the  wavelet  coefficient  mag¬ 
nitudes,  3)  zeroing-out  all  coefficients  with  magnitudes  less  than  the  threshold  and 
retaining  the  others,  and  4)  inverse  transforming  the  thresholded  coefficient  set  to 
yield  the  denoised  time  domain  signal  [4,7,8,14,15,17-19,44,61].  The  effectiveness  of 
this  approach  is  based  on  selecting  a  threshold  value  that  1)  retains  coefficients  con¬ 
taining  a  majority  of  desired  signal  contributions  while  2)  zeroing-out  coefficients  that 
are  dominated  by  noise  contributions.  Due  to  the  compaction  property  of  the  wavelet 
transform,  there  are  relatively  few  large  magnitude  coefficients.  Thus,  a  majority  of 
the  coefficients  can  be  zeroed-out  which  minimizes  the  remaining  noise  contribution 
in  the  denoised  response. 

Summary 

This  chapter  presented  the  relevant  technical  background  information  on  Burst 
Detection  and  Location,  RF  Fingerprinting,  DT-CWT,  and  Denoising.  The  informa¬ 
tion  here  supports  subsequent  material  presented  in  the  document. 


25 


III.  Methodology 

This  chapter  provides  the  overall  demonstration  process  used  for  generating  results 
and  conducting  analysis.  A  detailed  description  is  included  for  both  the  “Signal 
Collection”  (Section  [372])  hardware  and  “Post-Collection  Processing”  (Section  13. 3p 
software  processes.  The  primary  hardware  used  for  signal  collection  was  AFIT’s  RF 
Signal  Intercept  and  Collection  System  (RFSICS)  with  subsequent  data  processing 
accomplished  exclusively  in  a  MATLAB®  environment.  Denoising  using  the  DT- 
CWT  is  described  in  Section  13.41  Threshold  determination  for  the  various  processes 
is  described  in  Section  13.51 

3.1  Overall  Demonstration  Process 

Figure  [3711  shows  the  overall  demonstration  process  that  was  used  for  generating 
all  results  presented  in  Chapter  4.  The  dashed  boundaries  denote  the  processes 
that  are  primarily  conducted  in  hardware  and  software.  The  “Signal  Collection” 
hardware  process  consisted  of  placing  communication  devices  the  RFSICS  in  a  given 
electromagnetic  environment  and  making  signal  collections.  The  collected  signal  data 
(a  series  of  complex  valued  samples)  is  passed  along  for  subsequent  Post-Collection 
Processing  which  was  accomplished  exclusively  in  a  MATLAB®  environment.  The 
implementation  and  functionality  of  various  processes  in  Figure  13.11  is  discussed  in 
the  following  sections. 

3.2  Signal  Collection  Process 

Classification  performance  was  demonstrated  for  two  cases,  including:  1)  Intra¬ 
manufacturer  where  all  devices  are  from  a  given  manufacturer  and  have  different  serial 
numbers,  and  2)  Inter-manufacturer  where  at  least  one  of  the  devices  is  from  a  dif¬ 
ferent  manufacturer.  A  summary  of  manufacturers,  device  serial  numbers  and  signals 
considered  is  provided  in  Table  13.11  Consistent  with  the  overall  research  objective, 
the  table  shows  that  results  were  not  generated  for  all  combinations  of  manufactur- 
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Figure  3.1:  Overall  demonstration  process  for  signal  collection,  analysis  signal  gen¬ 
eration,  burst  detection  and  start  location,  fingerprint  extraction,  and  classification. 


ers,  devices  and  signals.  Rather,  selected  combinations  were  used  to  generate  results 
to  sufficiently  support  final  research  conclusions-it  is  believed  that  results  from  an 
exhaustive  analysis  would  not  fundamentally  change  these  conclusions.  From  an  op¬ 
erational  perspective,  the  potential  number  of  combinations  that  may  be  of  interest 
to  the  broader  technical  community  is  nearly  limitless  and  based  on  tens  of  man¬ 
ufacturers,  tens  of  device  types  per  manufacturer,  and  hundreds  of  serial  numbers 
per  device  type.  Considering  various  combinations  of  alternatives  remains  an  area  of 
interest  for  future  research  and  is  subject  to  technical  community  interest. 

For  all  results  presented,  the  signals  were  collected  with  both  the  device  under 
test  and  the  RFSICS  in  an  anechoic  chamber.  Basic  functionality  of  the  RFSICS  is 
provided  by  Agilent’s  E3238S  system  [1],  This  includes  an  RF  front-end  collection 
range  of  20.0  MHz  to  6.0  GHz  from  which  a  band  of  interest  is  selected  using  a 
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Tabic  3.1:  Device  manufacturers,  serial  numbers,  and  signal  types  (802. 11 A  and 
802. 11G)  used  for  generating  Chapter  4  results. 


Manu 

Serial  Number  /  Signal  Type 

Cisco 

N4U9  /  A&G 

N4UD  /  A&G 

N4UW  /  A&G 

N4PX  /  A&G 

Linksys 

0306  /  A&G 

0307  /  A 

Net  gear 

0273  /  A 

0217  /  A 

Dell 

BTA4 / A 

Airmag 

2C01  /  G 

tunable  RF  filter  with  fixed  bandwidth  of  36.0  MHz.  The  selected  RF  band  is  down- 
converted  to  an  Intermediate  Frequency  (IF)  of  70.0  MHz  and  passed  to  a  digitizer. 
The  digitizing  process  consists  of  down-conversion  (near  baseband),  12-bit  analog- 
to-digital  conversion  at  95  M  samples-per-second  (sps),  digital  filtering  (user  defined 
bandwidth),  Nyquist  compliant  sub-sampling,  and  data  storage  as  complex  In-phase 
(I)  and  Quadrature  (Q)  components.  A  digital  filter  bandwidth  of  18.56  MHz  was 
selected  for  all  802.11A  signals  collected  for  this  work.  This  resulted  in  the  RFSICS 
automatically  applying  a  sub-sampling  factor  of  four,  for  a  final  sample  rate  of  fs  = 
23.75  Msps  and  corresponding  sample  interval  of  Ts  =  1  /  fs  ss  42.1  ?zsec  per  sample. 
The  typical  collected  SNR  for  the  chamber  collected  signals  is  on  the  order  of  SNR  = 

40  dB. 

3.3  Post- Collection  Processing 

Post-Collection  Processing  in  Figurel34~1is  accomplished  exclusively  in  a  MATLAB® 
environment  using  the  near-baseband,  complex  I-Q  data  from  RFSIC  collections. 
Post-collection  processing  includes  analysis  signal  generation,  burst  detection  and 
start  location,  statistical  fingerprint  generation  and  signal  classification.  The  func¬ 
tionality  and  implementation  of  each  of  these  processes  is  discussed  in  the  following 
subsections. 

3.3.1  Analysis  Signal  Generation.  The  first  post-collection  process  of  “Per¬ 
fect”  Burst  Extraction  uses  the  near-baseband,  complex  I-Q  data  from  the  RFSIC  col- 


lections.  Extraction  is  accomplished  through  a  combination  of  automated  amplitude 
threshold  detection  followed  by  visual  analysis  and  manual  alignment  to  accurately 
identify  the  sample  number  corresponding  to  the  burst  start.  The  extracted  burst 
responses  are  digitally  filtered  using  a  baseband  filter  and  power-normalized.  A  6th- 
order  Chebyshev  digital  filter  was  implemented  having  a  -3  dB  bandwidth  of  7.7  MHz. 
At  this  point,  the  sample  frequency  of  the  filtered  signal  is  fs  =  23.75  Msps  which  ef¬ 
fectively  represent  oversampling  by  a  factor  of  approximately  1.5  times  Nyquist.  Pro¬ 
vided  that  the  RFSICS  collection  and  subsequent  post-processing  is  identical  for  all 
signals,  it  is  reasonable  to  assume  that  “recording  coloration”  (variation  in  amplitude, 
phase  and/or  frequency  characteristics)  induced  by  the  RFSICS  and  post-processing 
prior  to  burst  start  location,  statistical  fingerprint  generation  and  signal  classification 
is  approximately  identical.  This  is  important  in  the  overall  process  and  ensures  that 
final  results  are  based  on  as  received  signal  characteristics  and  features  versus  being 
unduly  influenced  by  signal-dependent  collection  and  post-processing  coloration. 

The  desired  “Analysis  Signal”  is  intended  to  simulate  varying  SNR  conditions 
that  typically  exist  in  an  operational  environments.  This  signal  is  generated  by  adding 
like-Eltered,  power-scaled  noise  to  the  digitally  filtered,  power- normalized  signal.  This 
is  done  by  generating  random  complex  AWGN  that  is  filtered  using  the  same  digital 
Elter  as  used  for  the  signal.  The  filtered  noise  signal  is  then  power-scaled  to  achieve 
the  desired  analysis  SNR  when  added  to  the  filtered  signal.  A  representative  instan¬ 
taneous  amplitude  response  from  a  collected  802. 11 A  RF  burst  is  shown  in  Figure  [3721 
for  analysis  SNRs  of  SNR  =  10  dB  and  SNR  =  0  dB. 

3.3.2  Burst  Detection  and  Start  Location.  The  sequential  burst  detection 
and  burst  start  location  process  is  implemented  relative  to  what  may  occur  in  an 
operational  collection  system,  i.e.,  a  real-time  system  samples  the  environment,  de¬ 
tects  the  “presence”  of  bursts  and  locates  the  burst  start  point  (sample  number) 
within  the  turn-on  region.  This  process  was  functionally  implemented  in  the  Locate 
Burst  Start  block  in  Figure  [37TJ  The  specific  burst  detection  and  location  techniques 
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Time  (irsec) 


Figure  3.2:  Instantaneous  amplitude  responses  for  collected  802. 11A  signal:  (Top) 
Collected  Signal,  (Middle)  Filtered  Signal-plus-AWGN  at  SNR  =  10  dB  and  (Bot¬ 
tom)  Filtered  Signal-plus-AWGN  at  SNR  =  0  dB. 

that  were  implemented  in  this  block  include  Fractal-Bayesian  Step  Change  Detector 
(Fractal-BSCD)  ( Section [2. 1.11).  Traditional  Variance  Trajectory  (VT)  ( Section [2. 1.2p 
and  Denoised  VT  (Section [3. 41).  As  presented  in  Chapter  4,  results  were  generated  us¬ 
ing  each  of  these  techniques  to  characterize  1)  their  burst  detection  and  location  error 
performance,  2)  their  performance  relative  to  each  other,  and  3)  their  corresponding 
error  impact  on  subsequent  fingerprint  classification  performance.  It  became  evident 
throughout  the  research  that  reliable  comparison  of  error  impact  on  classification 
performance  could  only  be  accomplished  if  all  the  same  bursts  were  used  for  classifi¬ 
cation  following  detection  and  location.  To  ensure  a  fair  comparison,  the  concept  of 
“dual-convergent”  bursts  was  developed  as  explained  next. 

Undetected  bursts  are  those  which  are  actually  received  yet  their  presence  is  not 
declared.  Detected  bursts  are  those  which  are  received  and  their  presence  is  declared. 
The  focus  of  this  work  is  on  detected  bursts  with  subsequent  algorithmic  process¬ 
ing  used  to  determine  burst  start  location.  For  those  cases  where  the  burst  start 


30 


location  algorithm  does  not  converge  in  accordance  with  prescribed  criteria  (num¬ 
ber  of  iterations,  parametric  tolerance,  etc.),  the  detected  bursts  are  designated  as 
“non-convergent”  and  a  default  burst  start  location  value  assigned.  When  algorithm 
convergence  occurs,  the  bursts  are  designated  as  “convergent”  and  the  estimated  lo¬ 
cation  assigned.  When  algorithm  convergence  occurs  for  identical  bursts  with  two 
different  burst  location  techniques,  the  bursts  are  designated  as  “dual  convergent.” 

3.3.2. 1  Burst  Detection.  This  process  is  similar  to  coarse  burst  de¬ 
tection  that  is  accomplished  in  an  the  RF  environment  to  detect  the  presence  of  RF 
bursts.  The  input  analysis  signal  is  first  segmented  into  contiguous,  non-overlapping 
sub-sections  or  windows  such  that  Ns  =  Nw  in  (12.4)1 .  While  not  a  requirement,  non¬ 
overlapping  windows  are  used  to  minimize  processing  time.  This  has  the  disadvantage 
of  producing  coarser  estimates  of  where  the  actual  burst  response  starts,  while  at  the 
same  time  capturing  more  signal  power  within  each  window  and  improving  detectabil¬ 
ity.  For  all  results  presented  in  Chapter  4,  a  window  size  of  Nw  =  512  signal  samples 
(21.6  psec)  is  used. 

Two  burst  detection  methods,  Traditional  VT  and  Denoised  VT  as  described  in 
Section  12.1 .21  and  Section lT4lrespectivelv.  are  applied  to  the  windowed  signal  data  and 
an  a-priori  coarse  detection  threshold  tDet  used  to  declare  detection.  Once  a  coarse 
detection  occurs,  the  corresponding  segment  of  windowed  signal  data  is  passed  on  for 
start  location  determination  where  it  is  assumed  that  an  actual  burst  start  occurs 
within  the  window.  However,  as  with  all  coarse  signal  detection  approaches,  false 
alarms  can  occur  with  bursts  falsely  declared  present.  Coarse  detection  performance 
results  are  provided  in  Section  14.2.11 

3. 3. 2. 2  Burst  Start  Location.  This  process  is  similar  to  coarse  burst 
detection  in  that  the  Traditional  VT  and  Denoised  VT  techniques  are  reapplied  to 
determine  the  final  start  location.  In  addition,  the  Fractal-BSCD  technique  in  Sec¬ 
tion  [2TJ] is  considered  as  well.  For  the  Traditional  VT  and  Denoised  VT  techniques, 
the  precise  start  location  is  indicated  by  the  time  (sample  number)  at  which  an 
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abrupt  change  occurs  in  the  VTa  response  of  (12.31) .  For  the  Fractal-BSCD  technique 
the  precise  start  time  location  corresponds  to  the  time  (sample  number)  at  which 
a  maximum  occurs  in  the  a-posteriori  PDF  of  (12. 2H.  The  effectiveness  of  these  ap¬ 
proaches  is  based  on  an  implicit  assumption  that  bursts  of  OFDM-based  signals  can 
be  modeled  as  having  a  step  change  response  in/near  the  turn-on  transient  region. 
This  assumed  response  is  consistent  with  802. 11A  specifications  [29]  and  has  been 
successfully  exploited  in  related  research  [42,63,67]. 

For  all  three  techniques,  the  segment  of  windowed  data  that  is  passed  from  the 
coarse  detection  process  is  further  sub-segmented  using  much  narrower  and  highly 
overlapped  windows.  The  overlapping  windows  allow  for  better  location  accuracy  at 
the  expense  of  increased  processing  time.  For  this  work,  a  window  size  of  Nw  =  20 
samples  (0.84  psec)  is  used  with  a  shift  of  Ns  =  2  samples  (84.2  nsec )  between 
consecutive  windows. 

For  demonstrating  performance  of  the  Traditional  VT  and  Denoised  VT  tech¬ 
niques,  an  a-priori  location  threshold  (t jj0C  %  ioet)  is  used  to  automatically 
estimate  the  burst  start  location  based  on  a  significant  peak  response  occurring  in 
VTa  of  (12.31).  When  a  significant  peak  is  located  the  signal  is  passed  on  for  sub¬ 
sequent  fingerprint  generation.  In  some  cases  no  significant  peak  is  found  and  the 
algorithm  does  not  converge  to  a  solution.  This  non-convergent  condition  can  occur 
if  there  is  no  burst  present  (coarse  burst  detection  false  alarm)  or  if  the  threshold  is 
set  too  high  for  the  burst  under  evaluation.  There  are  two  options  for  dealing  with 
non-convergent  bursts,  including:  1)  the  burst  can  be  discarded  without  subsequent 
processing,  or  2)  a  default  start  location  value  can  be  assigned  and  subsequent  pro¬ 
cessing  performed.  In  an  operational  environment  where  the  system  has  access  to  a 
large  number  of  bursts,  discarding  burst  may  be  a  reasonable  choice  and  have  minimal 
impact  on  final  system  performance.  For  this  work  the  probability  of  coarse  detec¬ 
tion  is  effectively  100%  given  that  collected  signals  are  first  passed  through  “Perfect” 
Burst  Extraction  (via  a  visual  and  manual  inspection  of  each  burst)  according  to  Fig¬ 
ure  13.11  Given  this  and  data  collection  limitations,  a  default  location  is  assigned  to 
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non-convergent  bursts  that  produce  no  significant  peak  in  the  VTa{i)  response.  The 
default  location  time  (sample  number)  is  chosen  to  correspond  with  the  last  sample 
in  the  window  of  data  passed  by  the  coarse  detection  process.  In  presenting  results 
in  Chapter  4,  non-convergent  pulses  are  only  included  when  characterizing  detection 
and  start  location  error  performance  of  the  three  techniques  considered.  As  explained 
earlier,  they  are  not  included  when  assessing  the  impact  of  this  error  on  end-to-end 
signal  classification  performance.  A  performance  comparison  of  Traditional  VT  and 
Denoised  VT  burst  start  location  performance  is  presented  in  Section  14.2.2.41 

In  assessing  performance  of  the  Fractal-BSCD  technique  in  Section  [2.1.11  it  was 
found  that  there  were  no  non-convergent  bursts.  The  100%  convergence  of  Fractal- 
BSCD  processing  is  ensured  given  the  maximum  operator  in  (12. 2b.  Relative  to  the 
Traditional  VT  and  Denoised  VT  techniques,  this  could  be  an  operational  disadvan¬ 
tage  as  there  is  no  inherent  back-up  capability  for  detecting  bad  pulses  (false  alarms). 
A  performance  comparison  of  Fractal-BSCD  and  Traditional  VT  burst  start  location 
performance  is  presented  in  Section  14.2.21 

3.3.3  Statistical  Fingerprint  Generation.  Following  burst  detection  and 
start  location,  the  RF  statistical  fingerprints  are  generated  using  the  process  shown 
in  Figure  [373]  As  indicated  within  the  dashed  lines,  the  Characteristics  and  Statistics 
generating  functions  are  identical  for  both  the  time  domain  (TD)  and  wavelet  domain 
(WD)  techniques.  A  signal  region  of  interest  is  selected  from  the  input  analysis  signal 
and  parsed  into  a  predefined  number  of  sub-regions  for  fingerprint  generation.  For 
the  802.11A/G  signals  considered  here,  the  burst  preamble  is  the  region  of  interest. 
This  choice  was  based  on  1)  previous  works  which  successfully  exploited  the  preamble 
[20,54,55],  and  2)  the  preamble  sequences  being  identical  for  all  bursts  per  the  802.11 
standard  [29].  Figure  13.41  shows  the  modulated  signal  response  for  the  standard 
preamble  comprised  of  10  short  followed  by  2  long  symbols.  For  all  results  presented 
in  Chapter  4,  a  total  of  Nr  =  3  fingerprint  regions  were  used  as  highlighted  in 
Figure  l3~4l  The  three  different  fingerprint  regions  include  1)  the  first  8.0  psec  which 
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corresponds  to  ten  short  OFDM  symbols,  2)  the  last  8.0  /isec  which  corresponds  to 
two  long  OFDM  symbols,  and  3)  the  entire  16.0  /z sec  preamble  (both  short  and  long 
symbols) . 

For  TD  feature  classification,  the  centered  subregion  characteristics  are  calcu¬ 
lated  using  (12.60  (12.121)  and  statistical  classification  features  calculated  using  (I2.13p. 
(12.141).  and  (12. 15jl  for  each  resultant  characteristic  response.  The  resultant  TD  RF 
fingerprint  (feature  vector)  consists  of  27  total  features  per  collected  burst  (3  subre¬ 
gions  x  3  signal  characteristics  x  3  statistics).  The  TD  fingerprint  for  burst  b ,  from 
device  (class)  c,  in  subregion  r  is  given  by 

FrC  =  [  ^r(/)> 

7 r(a),  7 r(<j>),  7 r(f),  (3-1) 

ACy-(cz),  /tr(0),  Kr(y)  ]  j 

where  b  =  1,2,3,...,  iV^  with  iVj,  being  the  total  number  bursts,  r  =  1,2,3,..., 
with  Nr  being  the  total  number  of  subregions,  and  c  =  1,  2,  3  is  the  class  index. 
Considering  the  Nr  =  3  subregions  as  used  here,  the  composite  TD  classification 
feature  vector  (1  x  27)  is  formed  using  (13.11)  and  is  given  by 


F^’c  F?’c  Fg'0 


(3.2) 


For  WD  feature  classification,  the  processing  is  identical  to  TD  processing  ex¬ 
cept  that  a  Dual- Tree  Complex  Wavelet  Transform  (DT-CWT)  decomposition  is  per¬ 
formed  in  each  subregion.  As  depicted  in  Figure  I2T3T  the  DT-CWT  decomposes  each 
subregion  into  five  levels  associated  with  different  wavelet  scales.  The  “complex  WD 
signal”  samples  are  calculated  using  (12.21  jl ,  followed  by  characteristic  generation  and 
centering  using  (12 .6j)— (12. 12j) .  The  statistical  classification  features  are  calculated  using 
(12.131).  (12.141).  and  (j2.15j).  The  resultant  WD  RF  fingerprint  (feature  vector)  consists 
of  135  total  features  per  collected  burst  (3  subregions  x  5  DT-CWT  decomposition 
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Figure  3.3:  Generation  process  for  statistical  RF  fingerprints.  The  Characteristics 
and  Statistics  generating  functions  are  identical  for  both  the  TD  and  WD  techniques 
and  implemented  using  fl2.6p-fj2.12h  and  fl2.13p-fl2.15l).  respectively  [32], 


Entire  802.11a  Preamble 

◄ -  (  16  msec  Duration  )  - ► 

Fingerprint  Region  3 


Figure  3.4:  802. 11 A  preamble  structure  showing  OFDM  modulated  signal  response 

and  fingerprint  regions. 


35 


levels  per  subregion  x  3  signal  characteristics  x  3  statistics).  Paralleling  the  TD 
development,  the  WD  fingerprint  for  burst  b,  from  device  c,  in  subregion  r  which  has 
been  decomposed  into  l  DT-CWT  levels  is  given  by 


<Tr,z(a),  <7r,j(0),  <^,z(/), 

7r,l(a)j  'Yr,l((fi)i  7r,z(/)j 

«r,z(a),  Kr>l(0),  Kr>i(f )  ]  , 


(3.3) 


where  l  —  1,2,3, ...  ,Ni  with  A)  being  the  total  number  of  DT-CWT  decomposition 
levels  per  subregion.  Considering  Nr  —  3  subregions  with  N/  —  5  levels  as  used  here, 
the  composite  WD  classification  feature  vector  (1  x  135)  is  formed  using  (13.31)  and  is 
given  by 


F  b’c  - 

r  WD  ~ 


Fb,C  T7\0,C  T7lO,C  T7»0,C  TT'O,' 

1,1  *  1,2  *  1,3  *  1,4  *  1,5 

Fb,c  u \b,c  -[7 \b,c  -r-yb,c  -r?b,c 
2,1  -T  2,2  *  2,3  *  2,4  ^  2,5 


5,c 


p6,c  jp6,c 


?b,c 


3,1  '  3,2  '  3,3  '  3,4 


3,5 


(3.4) 


3.3.4  MDA/ML  Signal  Classification.  Signal  classification  is  performed 
using  the  Multiple  Discriminant  Analysis/Maximum  Likelihood  (MDA/ML)  process 
described  in  Section  [2. 2. 31  For  all  MDA/ML  classification  results  presented  in  Chap¬ 
ter  4,  a  total  of  Nb  =  2000  bursts  were  used  from  Nd  =  3  different  802.11A/G  devices, 
with  each  device  denoted  as  Class  A,  Class  B,  and  Class  C.  Fingerprints  from  each 
class  (device)  were  used  to  form  a  single  composite  fingerprint  matrix  for  classifica¬ 
tion.  As  indicated  in  the  following  expressions,  the  composite  matrix  is  formed  by 
vertically  concatenating  the  feature  vectors  for  either  TD  using  (13. 2 j)  or  WD  using 
(13. 4 j).  The  formation  of  these  matrices  can  be  represented  as 


Tpl,l 
*  TD 

-p2,l 

*  TD  ■ 

TT ,  1 

•  *  TD 

Tv  " 

F  td  — 

tD’2 

*  TD 

T?2’2 
*  TD  ■ 

TT'JVf,,2 

•  r  WD 

Tv 

ttU,3 
*  TD 

F2,3 

r  TD  . 

rf-Nj,, 3 

•  r  TD 

Tv 
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(3.6) 


F1,1 

r  WD 

F2,1 
r  WD  ' 

•  r  WD 

Tv  ' 

F  wd  = 

pi, 2 
r  WD 

F2’2 

r  WD  • 

pM,2 
•  r  WD 

Tv 

rpl,3 
r  WD 

T?2’3 

r  WD  ■ 

•  r  WD 

Tv 

where  Ty  is  used  here  to  denote  vector  transposition ,  i.e. ,  the  vectors  are  transposed 
with  the  order  of  elements  within  each  vector  maintained.  For  =  2000  bursts  per 
class,  the  resultant  composite  Ftd  matrix  has  dimension  6000  x  27  and  the  resultant 
composite  F wd  matrix  has  dimension  6000  x  135.  The  composite  fingerprint  matrices 
in  (13. 5j)  and  (13.61)  are  column-wise  (i.e.  per  feature)  centered  and  normalized  to  unit 
standard  deviation.  The  centering  and  normalizing  processes  only  aid  in  fingerprint 
visualization  and  do  not  impact  subsequent  MDA/ML  classification  performance. 
The  impact  of  feature  selection  (TD  and  WD)  on  signal  classification  performance  is 
demonstrated  using  the  resultant  centered  and  normalized  RF  fingerprints  input  to 
the  MDA/ML  process. 

Monte  Carlo  simulation  and  K-fold  cross  validation  processes  are  used  with 
MDA/ML  signal  classification.  Monte  Carlo  simulation  is  used  to  ensure  statistical 
significance  and  K-fold  cross  validation  is  used  to  generalize  the  prediction  error  to 
an  independent  data  set  [25].  While  the  required  value  of  K  can  vary  as  a  function 
of  data  “behavior,”  values  of  K  =  5  and  K  =  10  are  common  choices  for  cross 
validation  [25].  Using  K  =  5  with  A/  =  2000  bursts  (fingerprints)  per  device,  the 
input  fingerprints  are  partitioned  into  K  =  5  equal  subsets  (400  each),  with  K  —  1  =4 
subsets  (1600  fingerprints)  used  for  training  and  the  remaining  “hcld-out”  subset  (400 
fingerprints)  used  for  classification  [25]. 

The  overall  process  for  MDA/ML  classification  with  K-fold  cross  validation 
is  shown  in  Figure  13.51  Accounting  for  a  total  of  NMq  independent  Monte  Carlo 
noise  realizations,  the  process  for  generating  average  classification  results  includes 
the  following  steps.  Note  that  the  Fold  Iteration  Accumulator  in  Figure  [3751  is  cleared 
prior  to  the  start  of  this  process. 
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Figure  3.5:  MDA/ML  classification  process  with  K-fold  cross  validation  [71]. 


1.  Generating  the  analysis  signal  for  a  given  SNR  per  Section  [3.3.11 

2.  Performing  burst  detection  and  start  location  per  Section  13.3.21 

3.  Generating  statistical  fingerprints  per  Section  13.3.31  for  the  technique  under 
evaluation  (TD  or  WD) 

4.  Generating  projection  matrix  W  per  (12. 16p  using  K  —  1  =  4  subsets  (80% 
of  the  fingerprints)  from  each  device  for  training  and  ML  classifier  parameter 
calculation 
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5.  Transforming  the  “held-out”  subset  (20%  of  the  fingerprints)  from  each  device 
as  “unknown”  inputs  using  W  and  classifying  each  per  ML  criteria 

6.  Accumulating  the  current  fold  classification  results 

7.  Selecting  the  next  K  —  1=4  blocks  for  the  next  fold 

8.  Repeating  Step  4  -  Step  7  for  K  —  1=4  additional  folds 

9.  Repeating  Step  1  -  Step  8  a  total  of  Nmc  times  using  different  independent 
AWGN  realizations  for  each  iteration  (Fold  Iteration  Accumulator  not  cleared) 

10.  Averaging  Fold  Iteration  Accumulator  results  to  obtain  average  classification 
performance  (Accounting  for  all  factors,  the  final  average  is  based  on  a  total  of 
Nmc  x  Aft  x  3  independent  classification  decisions.) 

11.  Repeating  the  process  for  each  desired  analysis  SNR 


Representative  MDA-transformed  training  fingerprints  and  trained  decision  bound¬ 
aries  calculated  from  ML  distributions  are  shown  in  Figure  [3jJaJ  for  802. 11A  signals 
at  SNR  =  40  dB.  The  corresponding  projection  of  “unknown”  MDA-transformed 


fingerprints  are  shown  in  Figure  l3dfb)  overlayed  with  trained  decision  boundaries 
from  Figure l3T[l(a)|  Note  that  even  under  these  high  SNR  conditions  incorrect  classi¬ 
fication  is  possible.  For  example,  one  of  the  Class  C  (*  markers)  fingerprints  is  clearly 
projected  into  the  Class  A  (x  markers)  ML  decision  region  and  would  be  incorrectly 
classified. 


3-4  DT-CWT  Denoising  Process 

Denoising  is  accomplished  using  the  DT-CWT  described  in  Section  [2731  with  the 
process  illustrated  in  Figure  13171  The  complex  input  signal  f(n)  is  transformed  using 
the  DT-CWT  which  outputs  complex-valued  wavelet  coefficients  from  the  Treel  and 
Tree2  filter  banks.  These  outputs  are  combined  to  form  real- valued  coefficients  d(n) 
according  to 
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(a)  MDA/ML  Training:  Decision  Boundaries  Calculated  From  ML 
Distributions. 


x  Class  A  Point 


+  Class  B  Point 
*  Class  C  Point 


(b)  MDA/ML  Classification:  Projected  Fingerprints. 


Figure  3.6:  MDA/ML 
SNR  =  40  dB.  Lower  surface  of 
decision  boundaries. 


a)|  Training  and  (b)  Classification  for  802.11A  signals  at 
shows  MDA  fingerprint  projections  and  trained 
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Figure  3.7:  Denoising  process  using  the  DT-CWT  in  Section [2731  [31]. 


d(n)  =  \/|2>eel(n)|2  +  |Tree2(n)|2  .  (3.7) 

The  d(n)  coefficients  in  (13. 7[)  are  compared  with  the  denoising  threshold  tDN 
and  a  punctured  set  of  coefficients  d'(n)  produced  by  setting  all  coefficients  below  toN 
to  zero  and  retaining  those  above  tnN,  he.,  W  where  d(n')  <  t^N,  T reelin')  =  0  and 
Tree2{n')  =  0.  An  Inverse  DT-CWT  (IDT-CWT)  is  then  applied  to  d'{n)  to  produce 
the  denoised  complex  output  signal  g(n).  The  denoised  coefficients  are  subsequently 
processed  using  the  Traditional  VT  technique  in  Section  12.1.21  to  generate  Denoised 
VT  results. 

The  impact  of  denoising  is  demonstrated  by  comparing  Traditional  VT  results 
with  Denoised  VT  results  in  Figure  13.81  Note  that  the  circled  region  highlights  the 
burst  start  location  at  t  —  0.  The  representative  amplitude  response  |/(n)|  is  from  an 
802. 11A  burst  at  SNR  =  40  dB  and  is  identical  in  both  figures.  As  a  side  note,  the 
16.0  fj, sec  preamble  response  is  clearly  evident  in  the  amplitude  response.  This  burst 
was  processed  along  with  an  SNR  =  0  dB  scaled  version  to  generate  the  VT(n)  results 
shown  for  each  technique.  As  indicated,  both  techniques  produce  nearly  identical 


41 


VT(n)  responses  at  SNR  =  40  dB  with  a  clear  distinct  peak  coinciding  with  the  burst 
start  time  at  t  —  0.  The  effect  of  denoising  is  most  evident  in  the  SNR  =  0  dB  results 
by  comparing  the  VT(n)  responses  just  prior  to  t  =  0,  the  actual  burst  start  time. 
The  t  <  0  region  of  collected  signals  only  contains  background  noise  contributions. 
Upon  close  inspection  of  the  SNR  =  0  dB  responses  for  t  <  0,  it  is  evident  that 
DT-CWT  denoising  has  effectively  reduced  the  background  noise  response.  While 
both  VT(n)  responses  at  SNR  =  0  dB  have  a  peak  near  t  =  0,  only  the  Denoised  VT 
response  has  the  desired  step  change  response  that  is  required  for  effective  threshold 
detection  and  burst  location. 

3.5  Threshold  Determination  Process 

Three  distinct  threshold  values  are  required,  including:  1)  tr>et  for  coarse  burst 
detection  per  Sectionl3.3.2.11  2)  for  burst  location  per  Section  13.3.2.21  and  3)  tr>N 
for  denoising  per  Section  13.41  All  SNR  dependent  threshold  values  were  determined 
a-priori  based  on  noise-only  analysis  using  100,000  AWGN  realizations.  The  random 
noise  realizations  were  generated,  filtered,  and  scaled  for  the  desired  analysis  SNR. 
For  determining  tDet  and  tLOC  thresholds  the  resultant  colored  noise  was  analyzed 
using  appropriate  window  parameters  for  a  given  technique.  In  determining  t for 
DT-CWT  denoising,  the  resultant  colored  noise  was  transformed  by  the  DT-CWT 
and  coefficients  retained  for  threshold  determination.  In  all  cases,  results  from  the 
100,000  noise-only  iterations  were  histogrammed  and  the  threshold  value  empirically 
chosen. 

In  selecting  a  tLOC  value  for  burst  location,  a  trade-off  is  made  between  the 
number  of  early  burst  location  estimates  and  the  number  of  non- convergent  solutions 
produced  by  the  algorithm.  The  final  :Uoc  values  were  selected  to  ensure  that  both  of 
these  conditions  are  present  and  observable  in  the  data.  When  comparing  Traditional 
VT  and  Denoised  VT  performance,  the  tiQc  value  is  further  constrained  to  provide  a 
similar  number  of  early  burst  location  (10%)  for  both  techniques  to  illicit  a  more  fair 
comparison. 
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Time  (nsec) 


(a)  Traditional  VT. 


Time  (iisecs) 
(b)  Denoised  VT. 


Figure  3.8:  Instantaneous  amplitude  response  of  an  802. 11 A  burst  and 
tional  VT  and  (b)  Denoised  VT  responses  for  SNR  =  40  db  and  SNR  = 
circled  region  highlights  the  burst  start  location  at  t  —  0.  [31]. 


(a)  Tradi- 
0  dB.  The 
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For  DT-CWT  denoising,  the  value  of  tDN  is  empirically  chosen  and  based  on 
the  histogram  bin  value  below  which  95%  of  the  noise-only  values  occur.  The  value 
of  tnet  is  chosen  using  conventional  noise-only  analysis  of  Probability  of  False  Alarm 
(Pfa)  and  Probability  of  Detection  (Pd)  as  represented  on  a  Receiver  Operating  Char¬ 
acteristic  (ROC)  curve.  Results  of  this  analysis  are  reported  in  Section  14.2.11 

Summary 

This  chapter  provided  implementation  details  for  Signal  Collection  and  Post- 
Collection  Processing,  Statistical  Fingerprint  Generation,  MDA/ML  Signal  Classifi¬ 
cation,  DT-CWT  Denoising  and  Threshold  Determination.  The  results  from  imple¬ 
menting  these  processes  are  provided  in  Chapter  4. 
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IV.  Results 


This  chapter  provides  modeling,  simulation  and  analysis  results  that  were  generated 
using  the  processes  detailed  in  Chapter  3.  The  research  involved  hundreds  of  sim¬ 
ulations,  requiring  hundreds  of  hours  of  processing  time  in  some  cases.  For  brevity 
and  to  ensure  succinctness,  only  a  subset  of  representative  results  are  presented  from 
selected  scenarios  to  fully  support  key  research  findings  and  contributions.  Results 
for  each  contribution  area  introduced  in  Section  11.21  are  presented  in  the  following 
subsections:  Bandwidth  Sensitivity  in  Section  14.11  Burst  Detection  and  Location 
in  Section  [4.21  MDA/ML  Classification  in  Section  H.31  and  Performance  Sensitivity 
Analysis  in  Section  14.41 

4-1  Bandwidth  Sensitivity 

Prior  to  assessing  burst  detection  and  device  classification  performance,  there 
was  one  important  parameter  that  needed  to  be  analyzed  -  the  post-collection  filter 
bandwidth  (BWpc).  As  shown  in  Figure  [3H]  of  Section  l3Tl  the  collected  burst  re¬ 
sponses  and  simulated  noise  are  digitally  filtered  prior  to  forming  the  desired  analysis 
signal.  In  previous  related  works  using  802.11  signals,  this  filter  bandwidth  was  simply 
fixed  at  a  reasonable  value  based  on  common  engineering  practice  [31,33,54,55]. 

Intra-manufacturer  classification  accuracy  using  three  Cisco  devices  is  presented 
versus  post-collection  filter  bandwidth  in  Figure  l4~T1  for  both  TD  and  WD  techniques 
using  802. 11A  signals  at  SNR  =  40  dB.  While  the  best  case  WD  classification  perfor¬ 
mance  is  approximately  2%  poorer  than  best  case  TD  performance,  the  WD  technique 
is  more  robust  and  classification  performance  varies  by  less  than  2%  over  the  range  of 
bandwidths  considered.  The  TD  technique  is  much  more  sensitive  and  exhibits  classi¬ 
fication  variation  of  nearly  6%,  with  poorest  TD  classification  performance  occurring 
at  BWpc  =  6.3  MHz. 

To  highlight  one  potential  cause  for  increased  TD  sensitivity,  a  few  filter  re¬ 
sponses  for  different  bandwidths  are  shown  overlayed  with  a  representative  802.11A 
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signal  PSD  in  Figure  14.21  The  three  filter  bandwidths  chosen  for  illustration  include 
BWpc  =  5.0  MHz,  BWpc  =  6.3  MHz  (worst  case  TD  classification  performance),  and 
BWpc  =  7.7  MHz.  Of  particular  note  is  how  each  of  the  filters  impact  the  802.11A 
OFDM  subcarrier  response  that  exists  near  7.5  MHz.  Clearly,  the  BWpc  =  7.7  MHz 
filter  effectively  passes  this  carrier  unaltered  while  each  of  the  other  two  filters  in¬ 
duce  some  degree  of  attenuation.  This  suggests  there  may  be  additional  informa¬ 
tion  in  the  higher  frequency  components  that  the  MDA/ML  classification  process 
is  more  effectively  exploiting.  However,  signal  attenuation  alone  cannot  account  for 
all  the  TD  performance  differences  in  Figure  l4~Tl  given  that  the  poorest  performing 
BWpc  =  6.3  MHz  filter  actually  attenuates  the  7.5  MHz  carrier  component  less  than 
the  better  performing  BWpc  =  5.0  MHz  filter  (-5.0  dB  versus  -16  dB).  Thus,  the 
filter  impact  on  noise  (attenuation  and  spectral  distribution)  must  be  considered  a 
contributing  factor  as  well. 

To  enable  comparison  of  both  techniques  at  their  best  performance  levels,  a 
post-collection  bandwidth  of  BWpc  =  7.7  MHz  was  used  for  generating  all  the  sub¬ 
sequent  burst  detection  and  device  classification  results  presented  in  Section  14.21  arid 
Section  14.31  respectively.  This  particular  bandwidth  choice  gives  the  TD  technique 
an  approximate  2%  advantage  in  device  classification.  This  will  be  considered  when 
presenting,  comparing  and  analyzing  subsequent  results. 

4-2  Burst  Detection  and  Location 

This  section  discusses  how  traditional  burst  detection  and  burst  start  location 
techniques  are  sensitive  to  varying  noise  conditions  and  how  this  sensitivity  impacts 
overall  classification  performance.  Analysis  indicates  that  improving  the  accuracy  of 
burst  detection  and  location  can  lead  to  improved  device  classification. 

4-2.1  Burst  Detection.  Receiver  Operating  Characteristic  (ROC)  curves 
were  generated  using  the  process  in  Section  13.3.2.11  to  characterize  performance  dif¬ 
ferences  between  the  two  burst  detection  techniques  -  Traditional  VT  and  Denoised 
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Figure  4.1:  Intra-manufacturer  classification  accuracy  versus  post-collection  filter 
bandwidth  for  TD  and  WD  techniques  using  802. 11A  signals  at  SNR  =  40  dB. 


Figure  4.2:  Overlay  of  representative  802. 11 A  signal  PSD  and  post-collection  filter 
responses  for  various  BWpc- 
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VT.  Results  in  Figure  [4731  show  that  at  SNR  =  6  dB  and  SNR  =  0  dB,  the  Denoised 
VT  technique  provides  a  higher  probability  of  detection  {Pd)  for  a  given  probability 
of  false  alarm  {Pfa)-  With  a  higher  Pd  for  a  given  Pfa,  the  Denoised  VT  technique 
detects  and  outputs  more  bursts  for  subsequent  processing  when  compared  with  the 
Traditional  VT  technique.  With  more  bursts  being  detected  and  forwarded,  it  is  pos¬ 
sible  to  correctly  classify  the  device  in  less  time  and  have  a  higher  confidence  in  the 
classification. 

4-2.2  Burst  Start  Location.  To  isolate  the  effects  of  burst  location  accu¬ 
racy  from  the  effects  of  burst  detection  error,  the  802. 11A  RF  bursts  were  manually 
detected  prior  to  burst  location  analysis.  Thus,  there  is  no  noise-only  data  input  to 
this  process  to  generate  false  alarms  and  Pd  =  100%.  All  histogram  results  in  this 
section  share  two  common  attributes,  including:  1)  the  correct  burst  locations  occur 
at  t  —  0  sec  and  2)  the  default  non-convergent  solutions  occur  at  t  =16  /rsec  (see 
Section [3.3.2.21  for  discussion  on  non-convergent  solutions). 

4-2.2. 1  Channel  Noise  Variability.  These  results  illustrate  the  effect 
of  channel  noise  variation  for  a  given  802. 11A  RF  burst  and  200  AWGN  realizations 
that  are  generated,  filtered,  scaled  and  added  to  achieve  the  desired  analysis  SNRs. 
Fractal-BSCD  and  Traditional  VT  estimation  results  are  shown  in  Figure  l4~4l 

At  higher  SNRs  the  two  methods  perform  similarly  as  the  noise  power  varies, 
with  primary  differences  beginning  at  SNR  =  9  dB.  Fractal-BSCD  degradation  is 
directly  attributed  to  the  a-posteriori  PDF  degradation,  as  calculated  per  (12. 2  p  and 
shown  in  Figure  [Til  The  strong  peak  response  in  the  PDF  diminishes  and  becomes 
more  uniformly  distributed  as  noise  power  increases.  Traditional  VT  degradation 
is  attributed  to,  and  affected  by,  threshold  selection  criterion.  For  the  non-optimum 
method  implemented  here,  the  threshold  criterion  is  not  always  satisfied  and  a  default 
start  value  is  assigned  -  a  missed  detection  or  non-convergent  solution.  This  is  shown 
in  Figure  l4~ip(b)|  as  a  peak  forming  at  t  =16  /zsec.  The  number  of  missed  detections 
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Figure  4.3:  Probability  of  False  Alarm  (Pfa)  versus  Probability  of  Detection  (Pd) 
ROC  curves  for  Traditional  VT  and  Denoised  VT  techniques  at 
and  p)]  SNR  =  0  dB.  [31], 


SNR  =  6  dB 
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(a)  Fractal-BSCD. 


(b)  Traditional  VT. 


Figure  4.4:  Impact  of  Channel  Noise  Variation  on  burst  location  error  using 


(a)  Fractal-BSCD  and  (b)  Traditional  VT.  Histogram  for  200  independently  gen¬ 


erated,  filtered  and  scaled  AWGN  realizations  with  a  given  802. 11 A  RF  burst  [33]. 


50 


at  lower  SNRs  can  be  reduced  by  changing  the  threshold.  However,  this  also  reduces 
estimation  accuracy  and  precision  at  higher  SNRs. 


4 -2. 2. 2  Burst-to-Burst  Variability.  These  results  illustrate  the  effect 
of  burst-to-burst  variation  using  a  given  AWGN  realization  that  is  generated,  filtered 
and  scaled  to  achieve  the  desired  analysis  SNRs.  Results  for  200  collected  802.11A 
bursts  with  Fractal-BSCD  and  Traditional  VT  estimation  are  shown  in  Figure  l4~5l 


As  with  the  channel  noise  impact,  the  two  methods  perform  similarly  at  higher 
SNRs.  Differences  arise  at  lower  SNRs,  with  the  Traditional  VT  method  degrading 
as  before  and  producing  missed  detections.  The  missed  detections  are  shown  in  Fig¬ 


ure  [O  Kb)  as  a  peak  forming  at  t  =16  /isec.  The  Fractal-BSCD  response  degrades 
differently  than  before,  becoming  multi-modal  at  lower  SNRs  and  producing  a  sig¬ 
nificant  number  of  detections  in  the  noise-only  portion  of  the  signal.  The  modes  are 
attributable  to  anomalous  spikes  in  a  specific  noise  realization.  This  is  consistent  with 
results  in  [23]  and  [68]  given  that  BSCD  processing  is  most  effective  when  non-gradual 
parameter  changes  occur.  At  lower  SNRs  the  amplitude  change  is  too  gradual  in  some 
bursts  for  the  BSCD  method  to  reliably  detect  them. 


4-2. 2. 3  Combined  Noise-Signal  Variability.  These  results  illustrate 
the  combined  effects  of  channel  noise  and  burst-to-burst  signal  variability.  In  this 
case,  200  AWGN  realizations  were  generated,  filtered  and  scaled  for  each  SNR  and 
added  to  each  of  the  200  collected  802. 11 A  bursts  -  a  total  of  40,000  unique  AWGN 
realizations  per  SNR.  Results  for  Fractal-BSCD  and  Traditional  VT  estimation  are 
shown  in  Figure  [4761 


In  this  combined  channel  noise  and  burst-to-burst  variability  case,  the  channel 
noise  effects  are  dominant.  This  is  evident  in  that  channel  noise  effect  results  in 
Figure  14.41  are  nearly  identical  to  the  combined  effects  results  Figure  14.61  including 


the  missed  detections  shown  in  Figure  EOIb)  as  a  peak  forming  at  £  =16  jusec. 
At  higher  SNRs  the  two  methods  perform  similarly  as  the  noise  power  varies,  with 
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(a)  Fractal-BSCD. 


(b)  Traditional  VT. 


Figure  4.5:  Impact  of  RF  Burst  Variation  on  burst  location  error  using  (a)  Fractal- 


BSCD  and  (b)  Traditional  VT.  Histogram  for  200  collected  802. 11A  RF  bursts  and 
one  generated,  filtered  and  scaled  AWGN  realization  [33] . 
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(a)  Fractal-BSCD. 


(b)  Traditional  VT. 


Figure  4.6:  Impact  of  Combined  Channel  Noise  and  RF  Burst  Variation  on  burst 
location  error  using  (a)  Fractal-BSCD  and  (b)  Traditional  VT.  Histogram  for  200 
independently  generated,  filtered  and  scaled  AWGN  realizations  and  200  collected 
802. 11 A  bursts  [33], 
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primary  differences  beginning  at  SNR  =  9  dB.  Fractal-BSCD  degradation  is  directly 
attributed  to  the  a-posteriori  PDF  degradation,  as  calculated  per  (12.2ft  and  shown  in 
Figure  EIU  The  strong  peak  response  in  the  PDF  diminishes  and  becomes  more  uni¬ 
formly  distributed  as  noise  power  increases.  Traditional  VT  degradation  is  attributed 
to,  and  affected  by,  threshold  selection  criterion. 

Relative  to  Fractal-Bayesian  Step  Change  Detector  (Fractal-BSCD)  technique, 
burst  detection  and  location  performance  was  best  using  a  Traditional  Variance  Tra¬ 
jectory  (Traditional  VT)  technique  which  provided  results  that  were  consistent  with 
perfect  burst  estimation  performance  at  higher  SNRs  (10  <  SNR  <  30  dB).  However, 
performance  for  both  techniques  diverged  at  lower  SNRs  (—3  <  SNR  <  10  dB)  [33]. 
This  shortfall  provided  an  impetus  for  subsequent  burst  detection  research  aimed  at 
improved  performance  at  low  SNRs  [31]. 

4-2. 2. 4  Combined  Noise-Signal  Variability:  Denoised  VT.  As  demon¬ 
strated  in  the  previous  sections,  burst  start  location  error  for  Fractal-BSCD  and 
Traditional  VT  becomes  symptomatic  at  SNR  <  9  dB  and  there  is  room  for  im¬ 
provement  for  the  lower  SNR  range.  In  accordance  with  Section  13.41  the  Denoised 
VT  process  consists  of  denoising  the  bursts  with  a  DT-CWT  prior  to  calculating  the 
Traditional  VT. 

Unlike  results  in  Section  14.2.2. II through  Section  14.2.2. 31  which  were  presented  as 
3-dimensional  histograms  in  Figure  14.41  through  Figure  l4~6l  the  discernable  differences 
in  results  of  this  section  were  not  readily  apparent  when  presented  as  3-dimensional 
histograms.  Thus,  the  results  in  this  section  are  presented  as  2-dimensional  his¬ 
tograms  for  a  given  subset  of  SNRs  considered.  The  results  in  Figure  14771  show  the 
improvement  achieved  at  SNR  =  6  dB  and  SNR  =  0  dB  when  denoising  is  employed. 

For  the  SNR  =  6  dB  results,  the  Denoised  VT  technique  outperforms  the  Tra¬ 
ditional  VT  technique  by  1)  correctly  locating  24%  more  of  the  burst  start  locations 
while  2)  experiencing  a  tighter  distribution  near  the  main  peak  response.  Similar 
improvement  is  demonstrated  for  the  SNR  =  0  dB  results.  While  both  techniques 
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Figure  4.7:  Probability  Distribution  Functions  (PDF)  for  burst  start  location  error 
using  Traditional  VT  and  Denoised  VT  at  (a)  SNR  =  6  dB  and  (b)  SNR  =  0  dB  [31]. 
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experience  a  main  peak  that  is  late,  the  Denoised  VT  technique  correctly  locates  3.7% 
more  of  the  burst  start  locations  while  also  exhibiting  a  tighter  distribution  near  the 
main  peak  response. 

4-2.3  Error  Impact  on  Device  Classification.  In  an  operational  implemen¬ 
tation,  only  those  bursts  causing  location  convergence  according  to  Section  13.3.2.21 
would  be  used  for  further  processing.  Therefore,  for  comparing  classification  perfor¬ 
mance  only  “dual  convergent”  bursts  per  Section  13711  are  used,  i.e. ,  only  the  bursts 
that  result  in  a  converged  location  solution  from  both  techniques  being  evaluated. 
All  other  bursts  that  resulted  in  a  converged  location  solution  from  only  one  of  the 
techniques  are  excluded  from  subsequent  classification.  This  approach  was  adopted 
based  on  early  results  which  showed  that  singly  convergent  bursts  unduly  biased  re¬ 
sults  in  favor  of  the  technique  yielding  the  most  converged  solutions.  The  distribution 
differences  (and  their  associated  fingerprints)  account  for  the  only  differences  between 
the  two  techniques  being  processed  by  the  classifier.  Classification  results  in  this  sec¬ 
tion  were  generated  using  a  mix  of  manufactured  devices,  including  two  from  Cisco 
(N4U9  as  Class  A  and  N4UW  as  Class  B)  and  one  from  Dell  (BTA4  as  Class  C).  Given 
the  two  Cisco  devices  are  very  close  in  serial  number  their  discrimination  inherently 
presents  the  greatest  classification  challenge. 

4-2.3. 1  Fractal-BSCD  and  Traditional  VT  Classification.  Figure  l4~8l 
shows  average  MDA-ML  classification  accuracy  with  the  effects  of  Perfect,  Fractal- 
BSCD  and  Traditional  VT  burst  detection  error  included.  In  this  case,  Perfect  results 
are  obtained  using  a  start  location  based  on  visual  inspection  of  each  collected  burst. 
To  determine  if  perfect  burst  location  provides  best  possible  MDA-ML  classification 
accuracy,  a  uniform  randomly  distributed  error  was  added  to  perfect  start  location  es¬ 
timates  and  results  generated  for  comparison.  As  shown,  the  Perfect  with  Random  Er¬ 
ror  results  are  consistent  with  Perfect  results  and  marginally  better/poorer  for  SNR 
below/above  approximately  14  dB,  respectively.  With  respect  to  the  burst  location 
estimation  error  impact  to  classification  performance,  the  Traditional  VT  technique 
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Figure  4.8:  Average  MDA-ML  classification  accuracy  with  Perfect,  Fractal-BSCD 
and  Traditional  VT  burst  detection  error  included.  [33]. 

was  consistent  with  Perfect  estimation  for  6  <  SNR  <  30  dB  but  under  performed  for 
—3  <  SNR  <  6  dB.  Traditional  VT  also  provided  considerable  improvement  when 
compared  with  the  Fractal-BSCD  technique  at  lower  SNRs  (—3  <  SNR  <  18  dB), 
i.e.,  for  a  given  classification  accuracy  in  the  range  of  50%-80%  the  required  SNR 
for  Traditional  VT  is  3-6  dB  lower  than  what  is  required  for  Fractal-BSCD. 

Classification  performance  is  commonly  illustrated  using  a  confusion  matrix 
that  shows  the  percentage  of  time  a  particular  input  class  is  estimated  as  one  of  the 
possible  classes,  with  the  diagonal  entries  representing  correct  classification.  Table  I4~T1 
shows  the  classification  confusion  matrix  for  perfect  burst  location  results  in  Figure l4~8l 
at  SNR  =  30  dB.  As  indicated  by  off-diagonal  entries,  the  greatest  confusion  exists 
in  intra-manufacturer  classification  with  Class  A  and  Class  B  inputs  being  mostly 
confused  with  each  other.  The  Class  B  input  is  errantly  classified  as  Class  C  a  small 
percentage  of  the  time  and  the  Class  C  input  experiences  no  confusion.  Collectively, 
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Tabic  4.1:  Classification  confusion  matrix  for  perfect  burst  location  results  in  Fig¬ 
ure  14.81  at  SNR  =  30  dB. 


Class  Estimate 

Input  Class 

A 

B 

c 

A 

89.5% 

10.5% 

0.0% 

B 

10.0% 

89.5% 

0.5% 

C 

0.0% 

0.0% 

100.0% 

these  results  illustrate  that  the  most  stressing  classification  challenge  is  posed  for 
intra- manufacturer  discrimination  (the  two  Cisco  devices). 

4-2. 3. 2  Traditional  VT  and  Denoised  VT  Classification.  To  assess  the 
impact  of  DT-CWT  denoising,  Denoised  VT  classification  results  were  generated  for 
comparison.  These  results  are  presented  in  Figure  14.91  which  shows  average  MDA- 
ML  classification  accuracy  with  the  effects  of  Perfect,  Traditional  VT,  and  Denoised 
VT  burst  detection  error  included.  As  before,  the  Perfect  results  provide  an  upper 
bound  on  achievable  performance.  As  indicated,  Traditional  VT  and  Denoised  VT 
performance  is  similar  for  SNR  >  6  dB  and  SNR  <  —2  dB.  For  —1  <  SNR  <  5  dB, 
the  Denoised  VT  technique  outperforms  the  Traditional  VT  technique  and  provides 
an  average  improvement  in  classification  accuracy  of  1.75%.  Relative  to  results  for 
perfect  burst  detection  and  location,  the  Denoised  VT  process  achieves  nearly  34%  of 
the  available  performance  improvement-when  used  with  MDA/ML  processing,  there 
is  little  more  to  be  gained  in  overall  classification  performance  by  improving  burst 
detection  and  location  accuracy. 

Confusion  matrix  results  for  the  SNR  =  3  dB  data  points  in  Figure  14.91  are 
shown  in  Table  14.21  Two  things  are  evident  when  comparing  Traditional  VT  and 
Denoised  VT  results,  including:  1)  minimal  difference  in  Class  A  and  Class  B  perfor¬ 
mance,  and  2)  greatest  improvement  occurring  in  correctly  classifying  Class  C  which 
exhibits  a  6%  increase.  These  results  are  consistent  with  what  is  expected  when 
considering  “What  level  of  improvement  is  achievable?”  Assuming  Perfect  results 
represent  an  upper  bound,  achievable  improvement  is  determined  by  comparing  di- 


Figure  4.9:  Average  MDA-ML  classification  accuracy  with  Perfect,  Traditional 

VT  and  Denoised  VT  burst  detection  error  included.  Results  obtained  for  “dual 
convergent”  802. 11A  bursts  from  a  mix  of  Cisco-Cisco-Dell  devices  [31]. 

agonal  entries  in  Table  14.21  for  Perfect  and  Traditional  VT  techniques.  For  Class  A 
and  Class  B  devices,  there  is  only  a  l%-2%  margin  for  improvement  in  correct  clas¬ 
sification.  However,  there  is  a  12%  margin  for  improvement  in  Class  C  classification. 
Thus,  the  Denoised  VT  performance  improvement  of  6%  for  Class  C  represents  50% 
of  the  possible  improvement. 


4-3  MDA/ML  Device  Classification 

As  concluded  in  Section  14.2.3. 21  and  highlighted  by  results  in  Figure  l4~9l  there  is 
minimal  additional  improvement  that  can  be  made  in  end-to-end  device  classification 
by  considering  alternative  burst  location  techniques.  The  reader  is  reminded  here  that 
the  focus  of  this  research  is  on  pro  of- of- concept  demonstration  without  optimization 
for  real-time  implementation.  Thus,  there  may  be  alternate  burst  detection  techniques 
that  are  more  computationally  efficient  than  those  considered  here.  However,  their 
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Tabic  4.2:  MDA/ML  classification  confusion  matrix  for  various  burst  detection 

methods  at  SNR  =  3  dB  [31]. 


Perfect 

Class  Estimate 

Input  Class 

A 

B 

C 

A 

68% 

21% 

11% 

B 

31% 

44% 

25% 

C 

14% 

17% 

69% 

Traditional  VT 

Class  Estimate 

Input  Class 

A 

B 

C 

A 

67% 

22% 

11% 

B 

31% 

42% 

27% 

C 

22% 

21% 

57% 

Denoised  VT 

Class  Estimate 

Input  Class 

A 

B 

C 

A 

67% 

21% 

12% 

B 

30% 

43% 

27% 

C 

18% 

19% 

63% 

application  to  the  RF  fingerprinting  process  detailed  in  Figure  I3TT1  of  Section  I3TT1  is 
beyond  the  scope  of  this  research  and  remains  an  area  of  future  research. 

Given  the  burst  detection  capability  detailed  in  Section  14.2.3.21  and  the  inher¬ 
ent  robustness  of  the  MDA/ML  classification  process  described  in  Section  f3.3.41  the 
research  emphasis  shifted  toward  improving  device  classification  by  considering  alter¬ 
nate  RF  fingerprint  features.  More  specifically,  the  DT-CWT  process  in  Section  12.31 
that  was  used  for  Denoised  VT  burst  detection,  was  next  used  for  generating  finger¬ 
prints  according  to  Section  13.3.31  The  incorporation  of  a  DT-CWT  prior  to  statistical 
feature  calculation  is  functionally  illustrated  in  the  RF  fingerprinting  process  depicted 
in  Figure  [3731  For  comparative  assessment  and  clarity  of  presentation  in  this  section, 
results  based  on  DT-CWT  fingerprints  are  referred  to  as  Wavelet  Domain  (WD)  re¬ 
sults  while  all  other  results,  including  all  those  presented  in  previous  sections,  are 
referred  to  as  Time  Domain  (TD)  results. 
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Various  combinations  of  device  manufacturers  (Cisco,  Netgear,  Linksys,  and 
AirMagnet)  and  signals  (802. 11A  and  802. 11G)  are  considered  for  demonstration 
with  specific  stressing  cases  considered  and  analyzed.  Using  three  Cisco  devices, 
classification  results  are  generated  and  analyzed  to  demonstrate  serial  number  dis¬ 
crimination.  This  is  the  most  stressing  case  considered  and  is  denoted  throughout 
as  “intra- manufacturer”  discrimination.  Using  a  combination  of  devices  from  various 
manufacturers,  classification  results  are  generated  and  analyzed  to  demonstrate  what 
is  denoted  as  “inter- manufacturer”  discrimination.  For  comparative  analysis,  results 
are  generated  using  TD  and  WD  fingerprints  generated  from  identical  collected  sig¬ 
nals  with  identical  Monte  Carlo  noise  realizations  that  are  appropriately  filtered  and 
scaled  to  achieve  desired  analysis  SNRs.  This  enables  a  one-to-one  comparison  of  TD 
and  WD  classification  results,  with  a  performance  “gain”  defined  as  the  difference 
in  required  SNR,  expressed  in  dB,  at  a  given  classification  accuracy  level.  This  is 
analytically  expressed  as  SNRwd  —  SNRtd  at  a  given  classification  performance. 
For  tracking  performance  improvement  and/or  degradation  throughout  this  section 
of  the  document,  the  performance  gain  at  an  80%  classification  accuracy  level  is  used 
per  “reasonable”  criteria  detailed  in  Section  11.2.11  and  is  shown  in  the  figures  as  a 
circled  region. 

4-3.1  Statistical  Fingerprint  Features.  The  ability  to  visualize  fingerprint 
features  can  be  insightful  for  both  feature  selection  and  performance  analysis.  Two 
important  properties  that  fingerprints  should  posses  to  increase  overall  classifica¬ 
tion  performance  are  uniqueness  and  temporal/spectral  stability.  Greater  fingerprint 
uniqueness  across  devices  provides  greater  separability  and  improved  classification 
performance.  Temporal  and  spectral  stability  of  fingerprint  features  is  also  impor¬ 
tant,  especially  for  the  MDA/ML  training  and  classification  process.  Ideally,  the 
statistical  fingerprint  features  used  for  MDA/ML  training  and  classification  do  not 
differ  significantly.  Given  the  signal  collection  conditions  used  for  this  research,  the 
temporal  and  spectral  stability  of  fingerprint  features  is  nearly  the  best  that  can  be 
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expected.  The  2000  bursts  used  for  all  of  the  results  presented  here  were  collected 
over  a  relatively  short  time  interval  (typically  less  than  0.5  sec)  and  in  an  anechoic 
chamber  void  of  multipath  and  channel  fading  effects.  The  uniqueness  of  fingerprint 
statistical  features  and  degree  of  temporal  stability  can  be  illustrated  using  what  are 
called  “Distinct  Native  Attributes”  (DNA)  in  RF  Fingerprint  DNA  plots. 

The  uniqueness  of  fingerprint  statistical  features  is  illustrated  in  Figure  14.101 
and  Figure  14.111  These  RF  DNA  plots  were  generated  by  randomly  selecting  250  col¬ 
lected  bursts  for  each  device,  scaling  them  to  achieve  SNR  =  20  dB,  and  averaging 
the  corresponding  statistical  fingerprints.  For  visual  clarity,  the  average  fingerprint 
features  are  normalized  within  each  segment  where  the  y-axis  segment  numbers  corre¬ 
spond  to  the  nine  statistical  measures  defined  in  (13.1  ft  and  (13.311.  The  number  of  DNA 
markers  per  segment  is  different  for  TD  and  WD  fingerprints.  For  TD  fingerprints, 
the  number  markers  is  a  function  of  the  number  of  signal  regions  used  for  fingerprint 
generation  as  expressed  in  (13. 2 j).  For  WD  fingerprints,  the  number  of  markers  is 
a  function  of  the  number  of  signal  regions  and  DT-CWT  levels  used  for  fingerprint 
generation  as  expressed  in  (I3.4[) .  The  RF  fingerprints  in  Figure  14.101  are  from  one  man¬ 
ufacturer  (Cisco)  and  typical  of  what  is  used  for  intra-manufacturer  discrimination. 
The  RF  fingerprints  in  Figure  14.111  are  from  three  different  manufacturers  (Cisco, 
Linksys  and  Netgear)  and  are  typical  of  what  is  used  for  inter-manufacturer  discrim¬ 
ination.  Two  conclusions  are  readily  apparent  by  analyzing  results  in  Figure  14.101 
and  Figure  14. Ill  including:  1)  relative  to  intra-manufacturer  fingerprint  features, 
the  inter-manufacturer  fingerprint  features  exhibit  greater  uniqueness  across  devices, 
and  2)  relative  to  TD  fingerprints,  the  WD  fingerprint  features  exhibit  greater  unique¬ 
ness  across  devices.  Subsequent  results  in  this  chapter  show  that  greater  uniqueness 
translates  to  better  overall  classification  performance. 

The  temporal  stability  of  fingerprint  features  is  demonstrated  in  Figure  14.121 
through  Figure  14.141  These  RF  DNA  plots  were  generated  by  randomly  selecting  25 
collected  bursts  for  each  device,  scaling  them  to  achieve  SNR  =  20  dB,  and  gen¬ 
erating  the  corresponding  fingerprint  for  each.  As  before,  the  fingerprint  features 


62 


Cisco  N4U9  Cisco  N4UD  Cisco  N4UW 

(a)  TD  Fingerprints. 


Cisco  N4U9  Cisco  N4UD  Cisco  N4UW 

(b)  WD  Fingerprints. 
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Figure  4.10:  Intra-manufacturer  average  RF  fingerprint  DNA  plots  showing |(a)|TD 
and  (b)  WD  fingerprints  based  on  250  randomly  selected  bursts  at  SNR  =  20  dB. 
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Cisco  N4U9  Netgear  0273  Linksys  0306 

(a)  TD  Fingerprints. 


Cisco  N4U9  Netgear  0273  Linksys  0306 

(b)  WD  Fingerprints. 
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Figure  4.11:  Inter-manufacturer  average  RF  fingerprint  DNA  plots  showing  |(a)|TD 
and  (b)  WD  fingerprints  based  on  250  randomly  selected  bursts  at  SNR  =  20  dB. 
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Ref  T2  T4  T6  T8  T10  T12  T14  T16  T18  T20  T22  T24 
(a)  TD  Fingerprints. 
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(b)  WD  Fingerprints. 


Figure  4.12:  Temporal  TD  Fingerprint  Stability:  (a)  TD  and  (a)  WD  Fingerprints 
for  25  randomly  selected  bursts  from  Cisco  N4U9  device  at  SNR  =  20  dB. 
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Ref  T2  T4  T6  T8  T10  T12  T14  T16  T18  T20  T22  T24 
(a)  TD  Fingerprints. 


-0.4 


-0.3 


0.2 


0.1 


Ref  T2  T4  T6  T8  T10  T12  T14  T16  T18  T20  T22  T24 

(b)  WD  Fingerprints. 


Figure  4.13:  Temporal  TD  Fingerprint  Stability:  (a)  TD  and  (b)  WD  Fingerprints 


for  25  randomly  selected  bursts  from  Linksys  0306  device  at  SNR  =  20  dB. 
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(a)  TD  Fingerprints. 
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(b)  WD  Fingerprints. 


Figure  4.14:  Temporal  TD  Fingerprint  Stability:  (a)  TD  and  (b)  WD  Fingerprints 


for  25  randomly  selected  bursts  from  Netgear  0273  device  at  SNR  =  20  dB. 
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are  normalized  within  each  segment  for  visual  clarity.  The  left- most  “Ref”  finger¬ 
prints  are  the  corresponding  average  reference  fingerprints  taken  from  Figure  14.101 
and  Figure  14.111  as  appropriate.  These  are  provided  for  comparison  with  the  ran¬ 
domly  selected  test  “T,:  fingerprints  which  are  presented  sequentially  with  increasing 
time  order.  Note  that  the  25  randomly  selected  test  bursts  are  different  than  the  250 
bursts  used  to  generate  the  average  reference  fingerprint.  The  results  in  Figure  14.121 
through  Figure  14.141  clearly  illustrate  a  dissimilar  degree  of  stability  among  the  fin¬ 
gerprint  features  being  used.  Note  that  the  effects  of  temporal  stability  on  the  overall 
MDA/ML  classification  is  outside  the  scope  of  this  work  and  is  reserved  for  future 
research. 

4-3.2  TD  vs.  WD  Performance:  802.11 A  Signals.  Intra-manufacturer  clas¬ 
sification  is  demonstrated  using  four  Cisco  devices  transmitting  an  802.11A  signal, 
with  results  presented  for  all  permutations  of  devices  as  shown  in  Table  14.31  Subse¬ 
quent  intra-manufacturer  discrimination  is  then  demonstrated  using  Permutation  #1 
which  presents  the  “most  stressing”  conditions  for  classification.  As  indicated  in 
Table  l4~3l  the  most  stressing  permutation  uses  three  Cisco  devices  having  serial  num¬ 
bers  that  differ  in  only  the  last  digit.  Thus,  it  is  assumed  that  these  devices  have 
been  manufactured  using  identical  components,  from  identical  lots,  with  identical 
processes,  under  identical  environmental  conditions.  Thus,  discriminating  between 
these  devices  presents  the  most  stressing  case  for  classification. 

Sensitivity  to  serial  number  variation  is  illustrated  in  Figure  14.151  which  shows 
intra-manufacturer  classification  results  for  all  four  permutations.  The  mean  across 

Table  4.3:  802. 11A  Cisco  intra-manufacturer  permutations. 


Serial  Number 

Perm 

N4U9 

N4UD 

N4UW 

N4PX 

1 

X 

X 

X 

2 

X 

X 

X 

3 

X 

X 

X 

4 

X 

X 

X 
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all  four  permutations  is  shown  by  the  filled  markers.  The  resuls  for  both  TD  and  WD 
techniques  show  that  Permutation  #1  and  Permutation  #3,  which  both  include  Cisco 
devices  with  serial  numbers  N4U9  and  N4UW,  present  the  most  stressing  cases  and 
yield  the  poorest  results  for  nearly  all  SNR  values  considered.  As  with  all  previous 
results,  Permutation  #1  is  the  most  stressing  case  at  80%  classification  accuracy. 

The  mean  classification  results  in  Figure  14.151  are  presented  again  in  Figure  14.161 
for  closer  inspection.  While  both  techniques  perform  similarly  at  SNR  >  25  dB, 
the  WD  fingerprinting  technique  outperforms  the  TD  technique  at  the  lower  SNRs. 
As  highlighted  in  the  circled  region,  the  WD  fingerprints  achieve  80%  classification 
accuracy  at  SNR  &  11  dB.  This  represents  a  gain  of  approximately  7  dB  with  respect 
to  equivalent  TD  fingerprinting  performance. 

Classification  confusion  matrices  are  presented  in  Table  l4~4l  for  Permutation  #1 
of  the  Cisco  devices  for  signals  at  SNR  =  11  dB.  As  indicated  in  the  lower  comparison 
matrix,  WD  fingerprinting  provides  improved  classification  performance  across  all 
three  classes,  with  the  greatest  improvement  of  28.1%  obtained  in  correctly  classifying 
Class  B.  One  common  result  with  both  fingerprinting  techniques  is  that  Class  A  and 
Class  C  devices  are  more  confused  with  each  other  and  confused  less  often  with 
Class  B.  With  respect  to  the  device  serial  numbers,  Class  A  and  Class  C  are  closer 
to  each  other  than  either  one  is  to  Class  B. 

Inter-manufacturer  classification  is  demonstrated  using  two  devices  each  from 
Cisco,  Netgear,  and  Linksys  transmitting  an  802. 11 A  signal,  with  results  presented 
for  device  permutations  shown  in  Table  l4~5l  Average  classification  performance  across 
all  device  permutations  are  shown  in  Figure  14.171  for  both  TD  and  WD  fingerprinting. 
While  both  techniques  perform  similarly  at  SNR  >  20  dB,  the  WD  fingerprinting 
technique  outperforms  the  TD  technique  at  the  lower  SNRs.  As  highlighted  in  the 
circled  region,  the  WD  fingerprints  achieve  80%  classification  accuracy  at  SNR  & 
2  dB.  This  represents  a  gain  of  approximately  5  dB  with  respect  to  equivalent  TD 
fingerprinting  performance. 
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Figure  4.15:  Intra-manufacturer  MDA/ML  classification:  Average  performance  for 
all  four  permutations  of  four  Cisco  devices  transmitting  802. 11A  signals. 


Figure  4.16:  Intra-manufacturer  MDA/ML  classification:  Average  performance 

across  four  permutations  of  four  Cisco  devices  transmitting  802. 11 A  signals. 
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Tabic  4.4:  Intra-manufacturer  confusion  matrices  for  TD  and  WD  fingerprinting: 
Permutation  #1  from  Table l4~3l with  802. 11A  signals  at  SNR  =  11  dB.  The  difference 
in  performance  between  the  two  techniques  is  provided  for  comparison. 


TD 

Class  Estimate 

Input  Class 

A 

B 

C 

A 

49.4% 

17.3% 

33.3% 

B 

18.5% 

65.9% 

15.6% 

C 

34.2% 

12.1% 

53.6% 

WD 

Class  Estimate 

Input  Class 

A 

B 

C 

A 

69.5% 

5.9% 

24.5% 

B 

5.3% 

94.0% 

0.7% 

C 

21.5% 

1.3% 

77.2% 

WD  -  TD 

Class  Estimate 

Input  Class 

A 

B 

C 

A 

20.1% 

-11.4% 

-8.8% 

B 

-13.2% 

28.1% 

-14.9% 

C 

-12.7% 

-10.8% 

23.6% 

Classification  confusion  matrices  are  presented  in  Table  [4761  for  Permutation  #1 
of  the  Cisco,  Netgear  and  Linksys  devices  at  SNR  =  2  dB.  Given  similar  results 
were  obtained  for  all  permutations  considered,  only  the  results  for  one  permutation 
are  presented  given  the  conclusions  drawn  are  generally  applicable  to  the  other  per¬ 
mutations.  While  the  WD  technique  increases  classification  performance  for  Cisco 
(Class  A)  and  Netgear  (Class  B)  devices,  there  is  a  decrease  in  Linksys  (Class  C) 
classification  performance  as  indicated  by  the  negative  diagonal  entry  in  the  lower  ma¬ 
trix.  The  greatest  improvement  of  30.2%  is  obtained  in  correctly  classifying  Class  A. 
Unlike  the  intra-manufacturer  discrimination  where  the  classes  are  similarly  confused 
regardless  of  the  fingerprint  technique,  the  inter-manufacturer  cross-class  confusion  is 
different.  The  TD  fingerprints  experienced  the  most  confusion  between  Class  A  and 
Class  B,  while  the  WD  fingerprints  showed  the  greatest  confusion  between  Class  B 
and  Class  C.  This  difference  accounts  for  the  greater  improvement  that  occurs  with 
Class  A. 
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Table  4.5:  802. 11 A  Inter-manufacturer  permutations. 


Cisco 

Netgear 

Linksys 

Perm 

N4U9 

N4UD 

0273 

0217 

0306 

0307 

1 

X 

X 

X 

2 

X 

X 

X 

Figure  4.17:  Inter-manufacturer  MDA/ML  classification:  Average  performance 

across  Cisco,  Netgear  and  Linksys  devices  transmitting  802. 11A  signals. 

4-3.3  TD  vs.  WD  Performance:  802.1 1G  Signals.  To  demonstrate  that  the 
classification  results  presented  up  to  this  point  are  not  unique  to  the  802. 11 A  signal, 
the  RF  fingerprinting  and  classification  process  was  applied  to  an  additional  OFDM- 
based  signal  to  demonstrate  broader  applicability.  This  was  easily  accomplished  using 
the  same  serial-numbered  devices  as  used  previously  by  operating  them  in  an  802. 11G 
signaling  mode. 

Using  the  same  four  Cisco  devices  (as  in  Section  14.3.21)  transmitting  an  802. 11G 
signal,  intra-manufacturer  discrimination  is  conducted  with  the  two  permutations 
shown  in  Table  l4~7l  Figure  14.181  shows  average  intra-manufacturer  classification  per¬ 
formance  across  the  two  permutations  of  Cisco  devices  for  TD  and  WD  fingerprinting. 
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Tabic  4.6:  Inter-manufacturer  confusion  matrices  for  WD  and  TD  fingerprinting: 
representative  permutation  of  devices  with  802. 11 A  signals  at  SNR  =  2  dB.  The 
difference  in  performance  between  the  two  techniques  is  provided  for  comparison. 


TD 

Class  Estimate 

Input  Class 

A 

B 

C 

A 

57.5% 

30.9% 

11.6% 

B 

34.7% 

53.5% 

11.7% 

C 

9.4% 

9.2% 

81.3% 

WD 

Class  Estimate 

Input  Class 

A 

B 

C 

A 

87.7% 

9.0% 

3.3% 

B 

7.5% 

71.5% 

20.9% 

C 

3.4% 

19.4% 

77.1% 

WD  -  TD 

Class  Estimate 

Input  Class 

A 

B 

C 

A 

30.2% 

-21.9% 

-8.3% 

B 

-27.2% 

18.0% 

9.2% 

C 

-6.0% 

10.2% 

-4.2% 

While  both  techniques  perform  similarly  at  SNR  >  20  dB,  the  WD  fingerprints  out¬ 
perform  the  TD  fingerprints  at  the  lower  SNRs.  As  highlighted  in  the  circled  region, 
the  WD  fingerprints  achieve  80%  classification  accuracy  at  SNR  «  11  dB.  This  rep¬ 
resents  a  gain  of  approximately  3  dB  with  respect  to  equivalent  TD  fingerprinting 
performance. 

Inter-manufacturer  classification  is  demonstrated  using  one  device  each  from 
Cisco,  Linksys,  and  AirMagnet  (shown  in  Table  14.81)  transmitting  an  802. 11G  signal 
Figure  14. 191  shows  average  classification  performance  for  TD  and  WD  fingerprinting. 
While  both  techniques  perform  similarly  at  SNR  >  20  dB,  the  WD  fingerprinting 
technique  outperforms  the  TD  technique  at  the  lower  SNRs.  As  highlighted  in  the 
circled  region,  the  WD  fingerprints  achieve  80%  classification  accuracy  at  SNR  ss 
2  dB.  This  represents  a  gain  of  approximately  2  dB  with  respect  to  equivalent  TD 
fingerprinting  performance. 
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Tabic  4.7:  802.1 


1G  Cisco  intra-manufacturer  permutations. 
Serial  Number 


Perm 

N4U9 

N4UD 

N4UW 

N4PX 

1 

X 

X 

X 

2 

X 

X 

X 

Figure  4.18:  Intra-manufacturer  MDA/ML  classification:  Average  performance 

across  two  permutations  of  four  Cisco  devices  transmitting  802. 11G  signals. 

4-3-4  Equivalent  TD  and  WD  Dimensionality.  Based  on  the  number  of 
classification  features,  the  WD  fingerprints  represent  an  approximate  5-fold  increase 
in  dimensionality  over  TD  fingerprints.  This  may  lead  one  to  conclude  that  the 
classification  improvement  with  WD  fingerprints  is  solely  attributable  to  using  an 
increased  number  of  features.  It  is  possible  that  the  performance  improvement  may 
be  the  result  of  more  exploitable  features  being  generated  from  the  DT-CWT  decom¬ 
position.  Thus,  it  is  reasonable  to  ask  “Is  the  noted  improvement  in  Section  14.3.21 
attributable  to  increased  feature  dimensionality,  more  exploitable  features,  or  both?” 
To  address  this  question,  results  were  generated  using  a  subset  of  27  selected  WD 
features  from  the  larger  135-feature  WD  fingerprints.  The  idea  was  to  compare  TD 
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Table  4.8:  802. 11G  Inter- manufacturer  permutations. 


Cisco 

Linksys 

AirMagnet 

Perm 

N4U9 

0306 

2C01 

1 

X 

X 

X 

Figure  4.19:  Inter-manufacturer  MDA/ML  classification:  Average  performance  us¬ 
ing  Cisco,  Linksys,  and  AirMagnet  devices  transmitting  802. 11G  signals. 


and  WD  performance  using  an  equivalent  number  of  features.  The  subset  of  WD 
features  was  selected  using  the  output  from  a  Generalized  Relevance  Learning  Vector 
Quantization  Improved  (GRLVQI)  classifier  [37-39].  The  GRLVQI  classifier  jointly 
selects  features  and  classifies  in  order  to  optimize  features  for  classification.  During 
this  process,  the  algorithm  calculates  and  outputs  a  relevance  rating  for  each  feature 
considered,  indicating  feature  importance. 

Using  WD  fingerprints  from  bursts  at  SNR  =  40  dB,  the  GRLVQI  classifier 
was  implemented  in  the  Waikato  Environment  for  Knowledge  Analysis  (WEKA)  en¬ 
vironment  [70]  and  used  to  determine  relevance  factors  for  all  135  WD  features.  The 
features  were  sorted  with  respect  to  their  relevance  and  the  27  most  relevant  features 
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Tabic  4.9:  Subset  of  27  most  relevant  WD  features  from  the  original  135  features. 
Relevance  ranking  (RR)  based  on  GRLVQI  classifier  output. 


RR 

Subregion 

WD  LVL 

Signal  Characteristic 

Statistic 

1 

Entire  Preamble 

4 

Amplitude 

Kurtosis 

2 

Short  Symbols 

4 

Amplitude 

Variance 

3 

Short  Symbols 

5 

Frequency 

Variance 

4 

Entire  Preamble 

4 

Amplitude 

Skewness 

5 

Short  Symbols 

1 

Amplitude 

Kurtosis 

6 

Entire  Preamble 

5 

Frequency 

Kurtosis 

7 

Short  Symbols 

3 

Amplitude 

Kurtosis 

8 

Long  Symbols 

2 

Phase 

Kurtosis 

9 

Entire  Preamble 

3 

Phase 

Kurtosis 

10 

Entire  Preamble 

3 

Phase 

Variance 

11 

Entire  Preamble 

1 

Frequency 

Variance 

12 

Short  Symbols 

3 

Amplitude 

Variance 

13 

Long  Symbols 

2 

Phase 

Skewness 

14 

Entire  Preamble 

5 

Amplitude 

Kurtosis 

15 

Entire  Preamble 

4 

Amplitude 

Variance 

16 

Entire  Preamble 

3 

Amplitude 

Kurtosis 

17 

Entire  Preamble 

4 

Frequency 

Kurtosis 

18 

Short  Symbols 

1 

Frequency 

Variance 

19 

Long  Symbols 

1 

Amplitude 

Kurtosis 

20 

Entire  Preamble 

5 

Phase 

Variance 

21 

Long  Symbols 

5 

Amplitude 

Variance 

22 

Short  Symbols 

2 

Amplitude 

Variance 

23 

Short  Symbols 

4 

Frequency 

Kurtosis 

24 

Entire  Preamble 

1 

Phase 

Variance 

25 

Entire  Preamble 

3 

Phase 

Variance 

26 

Long  Symbols 

1 

Phase 

Variance 

27 

Entire  Preamble 

1 

Phase 

Kurtosis 

retained  for  use  as  alternate  WD  fingerprints.  A  rank  ordered  listing  of  these  features 
is  provided  in  Table  14701  The  table  shows  the  final  relevance  ranking  (RR),  corre¬ 
sponding  preamble  subregion,  WD  level  (WD  LVL),  signal  characteristic  and  statistic 
for  each  ranked  feature,  ft  is  interesting  to  note  that  a  majority  of  the  most  relevant 
features  are  based  on  the  entire  preamble  region,  followed  by  the  variance  statistic 
and  then  a  tie  between  the  kurtosis  statistic  and  the  amplitude  characteristic. 
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Figure  4.20:  Inter-manufacturer  MDA/ML  classification:  Comparison  of  27-feature 
TD  and  27-feature  WD  performance  for  most  stressing  case  with  devices  transmitting 
802. 11 A  signals. 


The  27  most  relevant  WD  features  in  Table  [41)1  were  used  for  WD  fingerprinting 
and  performance  compared  with  27-feature  TD  fingerprinting  performance  under  the 
most  stressing  802. 11 A  intra-manufacturer  discrimination  case.  Figure  14.201  shows 
overall  classification  results.  As  highlighted  in  the  circled  region,  the  27-feature  WD 
fingerprints  achieve  80%  classification  accuracy  at  SNR  ps  19  dB.  This  represents  a 
gain  of  approximately  2  dB  with  respect  to  equivalent  27-feature  TD  fingerprinting 
performance.  Given  equal  dimensionality,  these  results  suggest  a  clear  increase  in 
exploitable  feature  information  using  the  DT-CWT  decomposition  process. 


4-4  Performance  Sensitivity  Analysis 

This  section  provides  results  that  address  classification  sensitivity.  Overall  ro¬ 
bustness  of  the  RF  fingerprinting  and  classification  process  is  assessed  for  three  spe¬ 
cific  cases,  including  variation  in  burst  location  error,  variation  in  MDA/ML  training 
and  classification  SNRs,  and  variation  in  MDA/ML  training  and  classification  signal 
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types.  Consistent  with  the  overall  proof-of-concept  research  objective,  the  results  here 
were  not  generated  with  a  goal  toward  achieving  optimal  performance.  Rather,  they 
address  a  few  of  the  most  apparent  “What  if?”  type  questions  that  are  of  interest  for 
operational  implementation  and  provide  a  basis  for  the  next  iteration  of  research. 

4-4--1  Effect  of  Burst  Location  Error.  The  effect  of  burst  location  error 
is  demonstrated  for  TD  and  WD  fingerprinting  using  random  burst  location  error. 
This  variation  addresses  the  operational  situation  where  equipment  used  for  collecting 
training  data  and  classification  data,  equipment  which  is  not  necessarily  co-located, 
may  be  operating  in  dissimilar  environments  that  are  less  than  ideal.  The  error 
considered  here  is  also  consistent  with  what  may  be  induced  by  laboratory  equipment, 
the  fidelity  of  which  can  impact  collected  signal  coloration  and  subsequent  burst 
location  accuracy.  Two  specific  random  error  distributions  are  considered,  including: 
1)  a  four-parameter  discrete  Beta  distribution  based  on  the  actual  observed  error  in 
post-processed  collected  data,  and  2)  a  uniform  distribution  having  minimum  and 
maximum  values  that  are  consistent  with  the  observed  error.  In  both  cases,  the 
location  error  is  randomly  applied  on  a  burst-by-burst  basis  to  the  perfect  burst 
location  data.  This  produces  what  is  referred  to  here  as  randomly  “jittered”  burst 
location  data. 

The  first  series  of  jittered  burst  results  was  generated  using  statistics  from  ob¬ 
served  location  error.  The  error  was  determined  on  a  burst-by-burst  basis  by  com¬ 
paring  sample  numbers  of  the  -3  dB  threshold  detected  bursts  and  the  corresponding 
manually  detected  perfect  bursts.  This  was  done  during  the  data  collection  and 
post-collection  processing  detailed  in  Section  13.21  The  resultant  histogram  for  ob¬ 
served  error  in  9134  collected  802. 11A  bursts  from  the  four  Cisco  devices  is  shown 
in  Figure  14.211  Based  on  statistics  of  the  observed  histogram  data  (mean,  standard 
deviation,  skewness,  and  kurtosis),  a  four-parameter  discrete  Beta  distribution  gen¬ 
erator  was  created  to  provide  simulated  location  error  similar  to  what  was  observed. 


78 


0.2 
0.18 
0.16 
0.14 

0.12 

LL 

Q  0.1 

CL 

0.08 
0.06 
0.04 
0.02 
0 

Figure  4.21:  Histogram  of  observed  burst  location  error  in  9134  collected  802. 11A 
bursts  from  four  Cisco  devices.  Simulated  error  results  for  the  four-parameter  discrete 
Beta  distribution  are  overlayed  for  comparison. 

Simulated  error  results  for  the  four-parameter  discrete  Beta  distribution  are  overlayed 
in  Figure  14.211  for  comparison. 

The  random  jitter  error  was  applied  to  perfect  burst  location  data  prior  to  ex¬ 
tracting  the  fingerprints  used  for  both  training  and  classification.  This  was  function¬ 
ally  implemented  within  Step  2  and  described  in  Section  13.3.41  Intra-manufacturer 
classification  results  (for  Permutation  ^2  in  Table  l4~3l)  using  observed  detection  error 
statistics  are  shown  in  Figure  14.221  for  both  WD  and  TD  fingerprinting  techniques. 
For  assessing  sensitivity  to  burst  location  jitter,  Figure  14.221  also  shows  performance 
for  perfect  burst  location  -  the  WD  technique  is  clearly  more  robust  than  the  TD 
technique.  Considering  the  circled  region  around  80%  classification  accuracy,  two 
conclusions  can  be  drawn:  1)  The  WD  technique  remains  superior  with  80%  classifi¬ 
cation  accuracy  achieved  at  SNR  sa  9  dB  for  both  jittered  and  perfect  burst  location 
error.  This  represents  gains  of  approximately  8  dB  (jittered)  and  6  dB  (perfect) 
with  respect  to  equivalent  TD  fingerprinting  performance;  2)  The  WD  technique  is 
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Location  Error  (samples) 
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less  sensitive  to  burst  location  error.  The  sensitivities  are  captured  by  considering 
the  SNR  differences  between  jittered  ( SNRj )  and  perfect  ( SNRp )  performance  at 
80%  classification  accuracy,  where  SNR&  =  SNRj  —  SNRp.  These  differences  are 
SNR&  ps  0  dB  for  WD  fingerprints  and  SNR&  ps  —2.0  dB  for  TD  fingerprints  where 
the  negative  sign  indicates  degradation.  The  near-zero  degradation  with  WD  finger¬ 
prints  clearly  indicates  the  WD  technique  is  more  robust  to  burst  location  error. 

The  second  series  of  jittered  burst  results  was  generated  using  uniformly  dis¬ 
tributed  error  of  ±6  samples  added  to  the  perfect  location  data.  This  particular 
range  of  values  was  chosen  based  on  the  maximum  observed  error  in  Figure  14.211  and 
presents  a  more  challenging  case  for  classification  (higher  mean  location  error  relative 
to  the  observed  statistics  case).  Results  in  Figure [47231  once  again  demonstrate  that 
WD  fingerprints  are  less  sensitive  to  location  error.  Comparison  of  WD  results  here 
with  those  in  Figure  14.221  shows  minimal  additional  degradation  with  uniformly  jit¬ 
tered  error.  However,  comparison  of  TD  results  here  with  those  in  Figure  14.221  shows 
considerably  more  degradation  with  uniformly  jittered  error. 

The  increased  sensitivity  is  captured  by  considering  the  SNR  difference  SNR& 
between  jittered  and  perfect  performance  at  80%  classification  accuracy.  The  differ¬ 
ence  for  WD  fingerprints  is  SNR&  «  —1  dB  which  is  marginally  different  from  the 
observed  jitter  case.  The  difference  for  TD  fingerprints  is  SNR&  ~  —12  dB  which  is 
twice  the  degradation  as  what  occurred  in  the  observed  jittered  case.  These  numbers 
indicate  that  WD  fingerprints  are  even  more  robust  than  previously  demonstrated 
with  observed  burst  location  error.  This  is  an  important  finding  for  two  reasons:  1)  it 
enables  subsequent  development,  demonstration  and  analysis  using  a  simple  uniform 
error  model  vice  requiring  a  rigorous  statistical  model  of  observed  location  error,  and 
2)  it  paves  the  way  for  subsequent  trade-off  studies  and  analysis  to  support  burst 
detector  selection  (hardware,  algorithm,  etc.)  for  system  implementation,  while  at 
the  same  time  addressing  the  question  “How  well  does  the  burst  detector  need  to 
perform?” 
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Figure  4.22:  Average  MDA-ML  classification  accuracy  for  802. 11 A  intra- 

manufacturer  discrimination  using  observed  burst  location  error  statistics. 


Figure  4.23:  Average  MDA-ML  classification  accuracy  for  802. 11 A  intra¬ 

manufacturer  discrimination  using  uniform  burst  location  error  statistics. 
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The  final  series  of  results  for  jittered  location  error  involves  the  use  of  dissimilar 
burst  location  accuracies  for  MDA/ML  training  and  classification  bursts.  The  intent 
is  to  represent  a  scenario  where  higher  fidelity  data  is  available  for  training  and  lower 
fidelity  data  is  used  for  classification.  The  assumption  is  that  higher  fidelity  data 
enables  better,  more  accurate  burst  location  while  lower  fidelity  data  yields  poorer, 
less  accurate  burst  location.  This  situation  may  occur  when  bursts  for  training  are 
collected  in  a  more  ideal  environment  and/or  with  better  equipment,  while  bursts  for 
classification  are  collected  under  poorer  environmental  conditions  and/or  with  poorer 
quality  equipment.  These  conditions  are  simulated  here  by  extracting  training  finger¬ 
prints  from  bursts  with  perfect  location  and  extracting  classification  fingerprints  from 
bursts  having  randomly  jittered  location  error.  In  this  case,  the  jittered  classification 
data  is  generated  using  the  statistical  distribution  of  the  observed  jitter  in  Figure  14.221 
with  a  variable  mean  delay. 

Classification  accuracy  for  intra- manufacturer  discrimination  is  shown  in  Fig¬ 
ure  14.241  for  a  mean  delay  of  0  to  90  samples  (0  to  3.79  /rsecs).  These  results  were 
generated  for  the  most  stressing  case,  Permutation  #1  in  Table  l4~3l  for  802. 11A  sig¬ 
nals  at  SNR  =  40  dB.  Note  that  performance  for  0  mean  delay  represents  an  upper 
bound.  As  indicated,  intra-manufacturer  discrimination  is  highly  sensitive  to  dissim¬ 
ilar  burst  start  location  error  with  performance  for  both  techniques  falling  below  80% 
accuracy  for  all  non-zero  mean  delay  values.  However,  the  WD  fingerprints  remain 
superior  for  a  majority  of  the  delay  values. 

Classification  accuracy  for  inter-manufacturer  discrimination  is  shown  in  Fig¬ 
ure  [4725]  for  a  mean  delay  of  0  to  90  samples  (0  to  3.79  /i secs).  These  results  were 
generated  for  Permutation  in  Table  14.51  for  802. 11A  signals  at  SNR  =  40  dB. 
As  indicated,  inter-manufacturer  discrimination  is  sensitive  to  dissimilar  burst  start 
location  error,  just  not  as  sensitive  as  intra-manufacturer  discrimination.  In  this  case, 
the  WD  fingerprint  performance  is  relatively  stable  for  mean  delays  below  14  samples 
(0.59  /xsecs)  while  the  TD  fingerprint  performance  immediately  decreases  over  this 
same  range.  Considering  the  circled  region  near  80%  classification  accuracy,  the  WD 


82 


Figure  4.24:  Average  MDA-ML  classification  accuracy  for  802. 11 A  intra- 

manufacturer  discrimination  using  dissimilar  burst  location  error. 


Figure  4.25:  Average  MDA-ML  classification  accuracy  for  802. 11 A  inter- 

manufacturer  discrimination  using  dissimilar  burst  location  error. 
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fingerprints  can  tolerate  up  to  55  samples  (2.32  /rsecs)  more  of  induced  mean  delay 
relative  to  the  TD  fingerprints. 

4-4-2  Effect  of  Dissimilar  Signal  SNRs.  A  comparison  is  made  here  between 
TD  and  WD  fingerprinting  performance  using  dissimilar  analysis  SNRs  for  MDA/ML 
training  and  classification.  Specifically,  the  training  burst  SNR  is  fixed  at  SNR  = 
40  dB  while  the  classification  burst  is  varied  at  SNR  <  40  dB.  Fingerprint  extraction 
and  classification  is  conducted  using  Permutation  ffl  in  Table  l4~3l 

Results  in  Figure  14.261  are  for  intra- manufacturer  classification  for  both  WD  and 
TD  fingerprints  using  dissimilar  analysis  SNRs  for  training  and  classification.  Rela¬ 
tive  to  performance  using  identical  training  and  classification  SNRs  (filled  markers), 
the  WD  technique  experiences  a  decrease  in  accuracy  for  all  SNR  <  30  dB  while 
the  TD  technique  actually  performs  better  at  SNR  >  18  dB  and  exhibits  decreased 
performance  at  SNR  <  18  dB.  However,  comparison  of  dissimilar  SNR  results  shows 
that  WD  performance  is  more  robust  for  SNR  <  20  dB.  As  highlighted  in  the  circled 
region,  WD  fingerprints  achieve  80%  classification  accuracy  at  SNR  ss  19  dB.  This 
represents  a  modest  gain  of  approximately  1  dB  with  respect  to  equivalent  TD  fin¬ 
gerprinting  performance.  This  is  approximately  7  dB  less  gain  when  compared  with 
performance  using  identical  SNRs  for  training  and  classification. 

Results  in  Figure  14.271  are  for  inter-manufacturer  classification  for  both  WD  and 
TD  fingerprints  using  dissimilar  analysis  SNRs  for  training  and  classification.  Unlike 
intra-manufacturer  results  which  exhibited  marginal  improvement  with  TD  finger¬ 
prints  over  a  limited  SNR  region,  there  is  only  degradation  in  the  inter-manufacturer 
results. 

Relative  to  performance  using  identical  training  and  classification  SNRs  (filled 
markers),  the  WD  technique  experiences  a  decrease  in  accuracy  for  all  SNR  <  12  dB 
while  the  TD  experiences  a  decrease  in  accuracy  for  all  SNR  <  20  dB.  Comparison 
of  dissimilar  SNR  results  shows  that  WD  performance  is  more  robust  overall  and 
performs  better  for  all  SNRs  considered.  As  highlighted  in  the  circled  region,  WD 
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Figure  4.26:  Average  MDA-ML  classification  accuracy  for  802. 11 A  intra¬ 

manufacturer  discrimination  using  dissimilar  SNRs. 


Figure  4.27:  Average  MDA-ML  classification  accuracy  for  802. 11A  inter¬ 

manufacturer  discrimination  using  dissimilar  SNRs. 
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fingerprints  achieve  80%  classification  accuracy  at  SNR  «  6  dB.  This  represents  a 
gain  of  approximately  6  dB  with  respect  to  equivalent  TD  fingerprinting  performance. 
This  is  approximately  2  dB  more  gain  when  compared  with  performance  using  iden¬ 
tical  SNRs  for  training  and  classification. 

4-4-3  Effect  of  Dissimilar  Signal  Types.  The  final  comparison  made  be¬ 
tween  TD  and  WD  fingerprinting  performance  involves  using  dissimilar  signal  types 
for  MDA/ML  training  and  classification  fingerprints.  Specifically,  training  finger¬ 
prints  are  generated  using  802. 11A  (802. 11G)  signals  with  classification  performed 
using  fingerprints  generated  from  802. 11G  (802. 11A)  signals.  Recall  that  the  col¬ 
lected  802. 11A  and  802. 11G  signals  are  from  the  same  physical  devices  operated  in 
two  different  modes.  Thus,  the  purpose  for  considering  dissimilar  signal  types  is  to  see 
if  there  are  inherent  signal  features  that  remain  unique  to  a  given  device  as  it  changes 
mode.  Fingerprint  extraction  and  classification  is  conducted  using  Permutation  ff T 
in  Table  14.31 

Results  in  Figure  14.281  are  for  intra- manufacturer  classification  for  both  WD 
and  TD  fingerprints  using  dissimilar  signal  types  for  training  and  classification.  For 
comparison,  classification  performance  is  shown  for  intra-manufacturer  discrimination 
of  802. 11A  signals  using  similar  signals  for  training  and  classification.  As  indicated 
by  the  encircled  data  points  at  SNR  =  40  dB,  the  intra-manufacturer  discrimina¬ 
tion  capability  is  very  poor  (50%  or  less)  using  either  WD  and  TD  fingerprinting 
techniques.  As  consistently  demonstrated  in  previous  sections,  the  WD  technique 
remains  more  robust  and  experiences  less  degradation  in  accuracy  when  compared  to 
the  TD  technique.  Given  these  intra-manufacturer  results  were  so  poor,  there  were 
no  additional  results  generated  for  inter-manufacturer  discrimination.  A  detailed  in¬ 
vestigation  into  the  cause (s)  of  such  poor  performance  was  not  within  the  scope  of 
this  research.  However,  there  are  two  issues  that  could  be  considered  a  good  starting 
point  for  such  an  investigation:  1)  The  same  hardware  devices  were  used  to  produce 
the  802. 11A  and  802. 11G  signals  which  fundamentally  operate  at  two  different  carrier 
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Figure  4.28:  Average  MDA-ML  classification  accuracy  for  intra-manufacturer  dis¬ 
crimination  using  dissimilar  signal  types  (802. 11A  and  802. 11G)  for  MDA/ML  train¬ 
ing  and  classification. 


frequencies.  Without  knowing  the  exact  device  details,  it  can  reasonably  be  assumed 
that  there  is  at  least  one  component  in  the  RF  transmission  chain  that  is  either  dif¬ 
ferent,  or  operated  differently,  between  the  two  modes  to  place  each  of  the  signals  at 
their  operating  frequencies.  Thus,  there  is  perhaps  dissimilar  coloration  that  impacts 
signal  features  such  that  they  are  not  the  same  across  the  two  operating  modes;  and 
2)  The  same  RFSICS  was  used  to  collect  the  two  signals.  Given  the  two  signals  are  at 
different  RF  carrier  frequencies,  the  internal  RFSICS  parameters  for  filtering,  down- 
conversion,  etc.,  are  necessarily  different  to  ensure  collected  signal  responses  reside  at 
baseband.  Thus,  there  is  perhaps  additional  coloration  due  to  RF /IF  collection  chain 
differences  in  the  RFSICS  that  can  further  impact  signal  features.  Collectively,  the 
RF  transmission  chain  of  the  devices  and  the  RF /IF  collection  chain  of  the  RFSICS 
could  be  inducing  unremovable  biases  in  collected  signals. 
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Summary 

This  chapter  provided  modeling,  simulation  and  analysis  results  that  were  gen¬ 
erated  using  the  processes  detailed  in  Chapter  3.  A  subset  of  representative  results 
were  presented  for  Bandwidth  Sensitivity,  Burst  Detection  and  Location,  MDA/ML 
Classification,  and  Performance  Sensitivity  Analysis.  Relative  to  corresponding  time- 
domain  (non-wavelet)  methods  and  results,  application  of  the  DT-CWT  provided 
improvement  for  all  burst  detection  and  RF  fingerprint  classification  scenarios. 


V.  Conclusion 


This  chapter  concludes  the  main  document  by  providing  an  overall  summary  of  re¬ 
search  activities,  a  summary  of  key  findings,  and  recommendations  for  subsequent 
research.  This  is  followed  by  an  appendix  that  provides  some  of  the  developmental 
MATLAB®  code  used  to  support  the  research. 

5.1  Research  Summary 

The  continued  proliferation  of  affordable  Radio  Frequency  (RF)  communica¬ 
tion  devices  has  greatly  increased  wireless  user  exposure  and  the  need  for  improved 
security  to  protect  against  spoofing.  Historically,  research  has  focused  on  the  detec¬ 
tion  and  mitigation  of  spoofing  using  bit-level  algorithmic  approaches.  More  recently, 
there  has  been  a  shift  toward  providing  added  security  within  the  Physical  (PHY) 
layer  of  the  Open  Systems  Interconnection  (OSI)  reference  model  by  exploiting  RF 
features  that  are  1)  inherently  unique  to  a  specific  device,  and  2)  are  difficult  to  repli¬ 
cate  by  an  unintended  party.  This  work  addresses  the  extraction  and  exploitation  of 
RF  “fingerprints”  to  classify  emissions  and  provide  hardware  specific,  serial  number 
identification-Specific  Emitter  Identification  (SEI).  The  related  SEI  concepts  that 
formed  the  foundation  for  this  research  are  collectively  embodied  in  previous  work 
on  RF  fingerprinting,  electromagnetic  signatures,  intrapulse  modulation,  and  unin¬ 
tentional  modulation  [11,23,24,30,34,51,64,66,68], 

Radar  systems  have  been  identified  using  SEI  techniques  that  exploit  inherent 
signal  features  that  are  unique  to  a  given  system  [9,40,60].  The  set  of  exploitable 
inherent  features  may  contain  unintentional  modulation  contributions  that  can  be 
influenced  by  any  number  of  environmental  and/or  hardware  issues,  some  of  which 
include  poor  system  design  (device  incompatibility),  improper  operation  ( over /under 
voltage),  and  physical  device  limitations  (operating  temperature  range)  [30,34,68]. 
Many  of  the  observed  unintentional  radar  modulation  effects  are  similar  to  what  exist 
in  modern  wireless  communication  systems  using  burst-like  waveforms.  This  begs 


the  question:  “Can  existing  SEI  methods  be  employed  with  wireless  communication 
signals  to  achieve  radar-like  SEI  capability?”  Answering  this  provided  the  motivation 
for  applying  a  Dual- Tree  Complex  Wavelet  Transform  (DT-CWT)  to  improve  burst 
detection  and  RF  fingerprint  classification. 

Despite  the  wealth  of  previous  work  that  forms  the  basis  for  this  research  [23, 
24,42,54,55,63,64,66,67],  the  task  of  automatically  detecting,  identifying  and  lo¬ 
cating  RF  communication  devices  remains  a  challenging  problem.  The  work  here 
addressed  four  main  aspects  of  this  problem,  including:  1)  the  selection  and  gen¬ 
eration  of  fundamental  signal  characteristics  (amplitude,  phase,  and/or  frequency), 
2)  the  feasibility  and  repeatability  of  detecting  and  locating  the  start  of  a  burst  using 
selected  waveform  feature(s)  amidst  channel  noise,  3)  the  identification  and  robust 
extraction  of  distinguishable  hngerprints-features  that  uniquely  characterize  the  un¬ 
intentional  modulation  of  a  device,  and  4)  the  performance  of  signal  classification 
under  varying  channel  conditions  and  Signal-to-Noise  Ratio  (SNR).  As  summarized 
below,  various  contributions  were  derived  from  the  research  while  addressing  each  of 
these  aspects. 

1.  SNR  Sensitivity  Analysis  [33]:  Except  for  the  two  most  relevant  earlier 
works  [54,55],  prior  works  lacked  a  detailed  sensitivity  analysis  of  burst  de¬ 
tection  and  fingerprint  classification  performance  under  varying  channel  SNR 
conditions.  To  address  this  deficiency,  this  work  analyzed  performance  of  two 
burst  detection  techniques,  including  the  Fractal-Bayesian  Step  Change  De¬ 
tector  (Fractal-BSCD)  and  Traditional  Variance  Trajectory  (Traditional  VT). 
With  respect  to  burst  location  estimation,  both  Fractal-BSCD  and  Traditional 
VT  techniques  provided  results  that  were  consistent  with  perfect  burst  location 
at  higher  SNRs  (10  <  SNR  <  30  dB).  However,  performance  for  both  tech¬ 
niques  diverged  at  lower  SNRs  (—3  <  SNR  <  10  dB).  With  respect  to  the  burst 
location  estimation  error  impact  to  classification  performance,  the  Traditional 
VT  technique  was  consistent  with  perfect  estimation  for  (6  <  SNR  <  30  dB) 
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but  under  performed  for  (—3  <  SNR  <  6  dB).  Traditional  VT  also  pro¬ 
vided  considerable  improvement  relative  to  the  Fractal-BSCD  at  lower  SNRs 
(—3  <  SNR  <  18  dB),  i.e.,  for  a  given  classification  accuracy  in  the  range  of 
50%-80%  the  required  SNR  for  Traditional  VT  is  3-6  dB  lower  than  what  is 
required  for  Fractal-BSCD. 

2.  Burst  Detection  at  Lower  SNR  [31]:  To  improve  burst  detection  and  lo¬ 
cation  capability  at  lower  SNRs,  the  DT-CWT  was  used  to  “denoise”  signals 
prior  to  applying  Traditional  VT  burst  detection.  Results  for  this  new  Denoised 
VT  technique  are  more  effective  and  provide  more  robust  burst  detection  and 
location  at  lower  SNRs  (—3  <  SNR  <  10  dB).  Relative  to  results  for  perfect 
burst  detection  and  location,  the  Denoised  VT  process  achieves  nearly  34%  of 
the  available  performance  improvement-when  used  with  Multiple  Discriminant 
Analysis/Maximum  Likelihood  (MDA/ML)  processing,  there  is  little  more  to 
be  gained  in  overall  classification  performance  by  improving  burst  detection  and 
location  accuracy. 

3.  TD  Fingerprint  Classification:  Given  demonstrated  improvements  in  burst 
detection  and  location,  the  research  emphasis  shifted  to  improving  upon  pre¬ 
vious  Time  Domain  (TD)  RF  fingerprinting  performance  using  the  newly  de¬ 
veloped  Wavelet  Domain  (WD)  RF  fingerprinting  technique.  To  assess  relative 
TD-WD  classification  performance,  it  was  necessary  to  replicate  previous  TD 
processing.  Given  that  all  previous  TD  work  was  based  on  a  fixed  post-collection 
bandwidth  BWpc  ~  9  MHz,  an  appropriate  value  based  on  sound  engineering 
practice  versus  best  or  optimal  criteria,  a  sensitivity  analysis  for  varying  BWpc 
was  conducted  to  determine  the  best  choice.  For  collected  802. 11A  signals  at 
SNR  =  40  dB,  this  analysis  indicated  that  TD  performance  was  very  sensitive 
and  exhibited  classification  variation  of  nearly  6%  for  5  MHz<  BWpc  <  9  MHz, 
with  best  case  near  100%  accuracy  realized  for  BWpc  =  7.7  MHz  and  poorest 
performance  realized  for  BWpc  =  6.3  MHz.  Thus,  BWpc  =  7.7  MHz  was  used 
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for  all  results  obtained  here  which  are  generally  better  than  previous  TD  results 
in  [54,55]  based  on  BWpc  ~  9  MHz. 

4.  WD  Fingerprint  Classification  [31-33]:  The  newly  developed  WD  RF  fin¬ 
gerprinting  technique  uses  coefficients  from  a  DT-CWT  decomposition  to  en¬ 
hance  statistical  fingerprint  features  and  improve  overall  device  classification 
performance.  Its  performance  was  demonstrated  in  four  stages: 

(a)  A  BWpc  sensitivity  analysis  was  conducted  similar  to  what  was  done  for 
TD  classification.  This  analysis  revealed  that  WD  performance  was  nearly 
insensitive  to  BWpc,  with  nearly  98%  accuracy  achieved  for  all  5  MHz  < 
BWpc  <  9  MHz.  For  comparative  TD-WD  assessment,  BWpc  =  7.7  MHz 
was  used  for  both  techniques  (best  case). 

(b)  Using  BWpc  =  7.7  MHz  (best  case  for  both  techniques),  improved  WD 
classification  performance  was  demonstrated  using  perfect  burst  location 
for  both  intra-manufacturer  (all  devices  from  the  same  manufacturer)  and 
inter-manufacturer  (a  mix  of  devices  from  different  manufacturers)  sce¬ 
narios  with  both  802. 11A  and  802. 11G  signals.  TD  and  WD  classification 
performance  was  compared  under  identical  scenarios  (devices,  signal  types, 
SNRs,  etc.)  using  a  “gain”  metric  defined  as  the  reduction  in  required  SNR 
for  the  WD  technique  to  achieve  the  same  classification  performance  as  the 
TD  technique.  For  80%  correct  classification  performance,  the  WD  tech¬ 
nique  provided  2  —  7  dB  gain  at  2  dB  <  SNRWD  <  11  dB.  The  approximate 
2%  best  case  TD  advantage  at  higher  SNRs  rapidly  diminishes  at  lower 
SNRs  that  are  more  consistent  with  operational  environments  [31,33]. 

(c)  The  previous  perfect  burst  location  results  were  based  on  27  TD  and 
135  WD  fingerprint  features.  To  address  the  question,  “Is  the  noted 
improvement  attributable  to  increased  feature  dimensionality,  more  ex¬ 
ploitable  features,  or  both?,”  results  were  generated  using  a  subset  of  27 
selected  WD  features  and  compared  with  27  feature  TD  results.  For  an 
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80%  classification  level,  the  WD  fingerprints  provided  a  gain  of  2  dB  at 
SNRwd  =  19  dB,  suggesting  a  clear  increase  in  exploitable  feature  infor¬ 
mation  in  the  DT-CWT  coefficients. 

(d)  Lastly,  WD  classification  performance  sensitivity  was  assessed  for  variation 
in  burst  location  error,  variation  in  MDA/ML  training  and  classification 
SNRs,  and  variation  in  MDA/ML  training  and  classification  signal  types. 
For  all  cases  considered,  the  WD  technique  proved  to  be  more  robust  and 
less  sensitive  when  compared  to  TD  technique. 

5.2  Recommendations  for  Future  Research 

As  noted  in  Section  11.1.21  the  choice  of  demonstrating  WD  fingerprinting  with 
OFDM-based  signals  was  motivated  by  two  factors,  including  1)  consistency  with 
previously  published  TD  work  [42,  54,  55,  63,  67]  and  2)  the  continued  emergence 
of  OFDM-based  signals  as  envisioned  for  4G  software  defined  and  cognitive  radio 
(SDR/CR)  communications  [21,26,48,72],  Relative  to  earlier  TD  work,  the  appli¬ 
cability  and  benefits  of  DT-CWT  fingerprint  features  has  been  clearly  demonstrated 
and  well-received  within  the  technical  community  [31-33].  However,  there  remains 
additional  topics  of  interest  that  could  be  investigated.  Some  of  the  most  evident 
include: 

1.  Optimization  of  Processes  or  Parameters:  As  used  for  demonstrating  DT- 
CWT  applicability  to  burst  detection  and  RF  fingerprinting,  there  are  numer¬ 
ous  processes  and  parameters  that  impact  performance.  Given  demonstration 
versus  optimization  was  the  goal  for  this  research,  the  degree  to  which  any 
given  factor,  parameter  and/or  combination  thereof  impacts  performance  was 
not  assessed.  As  developed,  implemented  and  demonstrated,  the  RF  finger¬ 
printing  process  is  well-suited  for  more  rigorous  optimization  using  a  Design  of 
Experiments  (DOE)  methodology  with  Analysis  of  Variance  (ANOVA).  The  op¬ 
timization  process  could  consider  any  number  of  conventional  techniques,  with 
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two  of  the  most  common  being  Genetic  Algorithms  (GA)  and  Response  Surface 
Methodology  (RSM). 

2.  Demonstration  Using  Different  Signals:  Demonstration  results  in  this 
research  were  based  on  collected  802. 11A  and  802. 11G  OFDM-based  signals. 
There  are  additional  OFDM-based  signals  that  are  emerging  for  4G  applications. 
For  example,  the  Worldwide  Interoperability  for  Microwave  Access  (WiMax) 
signal  has  emerged  and  is  rapidly  becoming  popular  for  establishing  “last  mile” 
communication  connectivity.  In  this  case,  the  designated  WiMax  base  station 
serves  a  similar  role  as  a  GSM  cellular  base  station  and  controls  user  activity 
within  a  defined  geographic  region.  Additional  work  could  be  performed  to 
address  the  use  of  RF  fingerprinting  to  provide  intra-cellular  WiMax  security. 

3.  Demonstration  of  Cross-Mode  Independence:  There  is  some  potential 
operational  benefit  if  “cross-mode”  device  discrimination  could  be  reliably  ac¬ 
complished,  i.e. ,  achieving  serial  number  SEI  based  on  fingerprint  features  that 
are  common  across  multiple  operating  modes  of  a  given  hardware  device.  For 
example,  there  are  IEEE  802.11  compliant  devices  that  support  multiple  modes 
(signal  types)  such  as  an  802.11A/B/G/N  device.  While  less  than  favorable,  the 
dissimilar  signal  type  results  in  Section  14. 4. 3lusing  802. 11A  and  802. 11G  suggest 
that  the  specific  features  considered  here  are  not  robust  enough  for  cross-mode 
classification.  Thus,  additional  cross-mode  work  could  be  performed  to  deter¬ 
mine  if  there  are  exploitable  underlying  RF  features  for  a  given  device  that  are 
independent  of  operating  mode. 

4.  Fused  Soft-Decision  Classification:  All  device  classification  results  pre¬ 
sented  in  this  work  are  based  on  averaging  what  may  be  called  “hard  decision” 
burst-by-burst  classification  decisions,  i.e.,  every  burst  input  to  the  MDA/ML 
process  is  associated  with  a  given  device  class  (Class  A,  Class  B,  or  Class  C) 
independent  of  how  other  input  bursts  are  classified.  In  many  applications  it  is 
often  possible  to  improve  performance  by  averaging  out  undesired  background 
noise  effects.  This  can  be  accomplished  in  various  system  processing  stages,  e.g., 
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RF,  IF,  pre-detection,  post-detection,  pre-classification,  post-classification,  etc. 
Given  the  communication  signals  of  interest  here  are  burst-like,  with  hundreds 
or  thousands  of  burst  generated  in  relatively  short  time  intervals,  it  is  reason¬ 
able  to  assume  that  device  classification  may  be  improved  using  what  may  be 
called  “soft-decision1'  device  classification.  Additional  work  could  consider  using 
knowledge  gained  by  analyzing  a  collection  of  multiple  burst-by-burst  classifi¬ 
cation  decisions  before  making  a  final  device  classification. 
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Appendix  A.  MATLAB ®  Code 

The  appendix  provides  the  main  MATLAB®  hies  used  to  functionally  implement  the 
processes  detailed  in  Chapter  3  and  used  for  obtaining  results  presented  in  Chapter  4. 
As  provided  below,  the  code  included  is  for  Burst  Detection  in  SectionfA.il  Preamble 
Location  in  Section  IA.21  Feature  Extraction  in  Section  IA.31  Device  Classification  in 
Section  IA.4I  arid  DT-CWT  Transformation  in  Section  [A~5l 

A.l  Burst  Detection 

Listing  A.l:  Code/Detect/PulseDetectV2.m 

7.  =  =  =  =  =  ====  =  =  =  =  =  =======  =  =  =  =  =  =  =====  =  =  =  ===  =  =  =  =  =  =  =======  =  =====  =  =  =  = 

•/.  Pulse  Detection  via  Amplitude  Thresholding 

7.  ============================================================= 

7. 

°l  Performs  Threshold  Amplitude  Detection  of  Pulses  and  Ouputs 
7,  a  Matrix  Containing  One  Pulse  Per  Row.  Amplitude  Detection 

7.  is  Accompli  shed  Using  a  Simple  Leading  Edge  Detector  Opera- 

70  tingon  a  Smoothed  Magnitude  Response  of  the  Input  Signal 

/  7. 

7.  function  [PlsMat  ,  PlsWdth  ,  PlsDb]  =  PulseDetectV2  (Z  ,  MaxPul  ,  .  .  . 

7.  AddS amp  ,  NumSmth  ,  PlsMin  ,  PlsMax  ,  Thresh  ,  NScr  ,  NPlot  ) 

7. 

7.  Created:  4  Nov  2008 
7.  By:  Dr.  Michael  A.  Temple 

/  Modified:  8  May  2009 
7.  By:  Dr.  Michael  A.  Temple 

7. 

/  Inputs7o 

7.  Z  =  Complex  Sampled  Input  Signal  (Column  or  Row  Vector)  .7. 

7.  MaxPul  =  Desired  Maximum  #  of  Pulses  to  be  Detected.  Actual 
7.  Number  in  Output  May  be  Less  Depending  on  the 

7.  the  Number  of  Detected  Pulses  Satisfying  PlsMin 

7o  and  PlsMax  Criteria7o 
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7,  AddSamp  =  Additional  #  of  Input  Samples  Included  Before 
7,  &  After  Threshold  Points  at  Edges  of  Pulse.0/, 

%  NumSmth  =  #  Samples  Smoothed/ Averaged  Across0/, 

7,  PlsMin  =  Min  #  Samples  in  Desired  Output  Pulse  Width 

7,  PI sMax  =  Max  #  Samples  in  Desired  Output  Pulse  Width0/, 

7,  Thresh  =  Desired  Detection  Threshold  Value  in  dB 
7,  Thresh  <  0  REQUIRED  !  °/„ 

7,  NScr  =  Output  Waitbar  Progress/Status  to  the  Screen? 

7.  1  =  Yes  0  =  No0/. 

7,  NPlot  =  Produce  Output  Plots? 

7,  1  =  Yes  0  =  No0/. 

7,  Guide  for  Selecting  Initial  Parameter  Values0/, 

7,  AddSamp:  Some  Number  <=  #  Samples  Between  Two  Closest  Spaced 
7,  Bursts  Divided  by  2. 

7,  NumSmth:  2°/,-5°/,  of  SHORTEST  Pulse  Duration.  Note  that  poorer 
7,  SNR  generally  requires  a  larger  NumSmth  value.0/, 

7. 

7,  Outputs 

7,  PlsMat  =  Output  Pulse  Matrix  with  One  Detected  Pulse  Per  Row. 

7,  When  Variable  Width  Pulses  are  Detected,  ALL 

7,  Non-MaxWidth  Pulses  are  Zero-Padded  in  last  Columns.0/, 

7,  PlsWdth  =  Pulse  Width  (#  Samples)  Between  Leading  &  Trailing 
7,  of  Detected  Pulse  Edges:  EstBetween  Leading  and 

7,  Trailing  Edges  of  the  Smoothed  Response.0/, 

7,  PlsDb  =  ACTUAL  Relative  Power  Level  (dB)  of  Output  Pulses  at 
7,  Leading  Edge  Detection  Point  of  SMOOTHED  Magnitude.0/, 

function  [PlsMat , PlsWdth , PlsDb]  =  PulseDetectV2 (Z , MaxPul , . . . 

AddSamp , NumSmth , PlsMin , PlsMax , Thresh , NScr , NPlot) 

LengthZ= length ( Z ) ; 

PlsMat  =  []  ; 

PlsWdth=  []  ; 

PlsDb=  []  ; 

7,  Ensure  /  Make  Z  a  Row  Vector 
Dim  =  size(Z); 
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if  Dim(2)==l  70  Column  Vector  Input  .  .  .  Change  to  Row  Vector 
Z  =  Z.J;  /  Use  Non-Conjugate  Transpose 

end 

7.  Pad/Extend  length  of  processed  ’TmpZ’  by  2*PlsMax  to  help 
%  mitigate  pulse  detection  issues  the  end  of  input  Signal  ’Z’ 
ExtZ=round (2* PlsMax ) ; 

TmpZ  =  [Z  ones ( 1 , Ext Z ) *min ( Z )] ; 

7.  Note:  Matlab’s  SMOOTH  Function  ALWAYS  returns  a  Column  Vector. 
7«  A  Transponse  is  Used  on  Smooth  Func  to  Restore  a  Row  Vector 
SmthZmag  =  20* loglO ( smooth ( abs ( TmpZ ), NumSmth) ’) ; 

TmpZmag  =  SmthZmag; 

LenTmpZmag  =  length ( TmpZmag) ; 

MinZ_Db=min  ( SmthZmag)  ;  °/0  Min  Value  of  Input  Signal 
MaxZ_Db=max ( SmthZmag) ;  7»  Max  Value  of  Input  Signal 
if  NScr==l 

BurstCons  =  waitbar (0 , ’ Starting  Burst  Detection  Loop’); 

end 

7«  Begin  Main  While  Loop 

7.  Initialize  Pulse  Detection  While  Loop  Variables 
PulseMatrix  =  [] ; 

PulseVec  =  []  ; 

MaxWidth  =  0; 

NumDet  =  0; 

PlsWdth=0 ; 

NPlsDet=0;  7»  Intialize  Pulses  Detection  Counter 
WhileMax  =  2*MaxPul  ;  7,  Set  Max  #  of  "While  Loop"  Iterations 
WhileCnt=0 ; 

while  NPlsDet  <  MaxPul  7»  Maximum  #  of  Pulses  to  be  Detected 
WhileCnt=WhileCnt+l ; 
if  WhileCnt  >=  WhileMax 
break 

end 

7.  Find  Smoothed  Peak  Response 
[MaxVal , MaxLoc ]  =  max (( TmpZmag) ) ; 


LowDex 


MaxLoc  ; 
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for  k 


=  l:2*PlsMax  %  Search  Left  to  Leading  Edge 
if(LowDex-l)  <  1  7«  First  Sample  Reached 
break  °/0  Stop  Search  ! 

else  °/0  Continue  Searching 

if  TmpZmag ( LowDex - 1 )  >  MaxVal  +  Thresh; 

LowDex=LowDex - 1 ;  7«  Index  #  at  Threshold 

else 

break ; 

end 

end 

end 

PlsLow  =  LowDex- AddSamp ; 
if  PlsLow  <  1 
PlsLow  = 1  ; 

end 

HghDex=MaxLoc ; 

for  k  =  l:2*PlsMax  7»  Search  Right  to  Trailing  Edge 
if(HghDex+l)  >  LenTmpZmag  7«  Last  Sample  Reached 
break  7,  Stop  Search  ! 

else  */.  Continue  Searching 

if  TmpZmag ( HghDex + 1 )  >  MaxVal  +  Thresh; 

HghDex  =  HghDex  +  l  ;  7«  Index  #  at  Threshold 

else 

break 

end 

end 

end 

PlsHgh  =  HghDex+AddSamp ; 
if  PlsHgh  >  LengthZ 

PlsHgh  =  LengthZ; 

end 

TmpWdth  =  HghDex  -  LowDex  ;  7«  Width  Between  Pulse  Edges 
7,  Check:  PlsMin  <  Temp  Width  <  PlsMax 
7.  Not  Satisfied  ->  Do  NOT  Include  Current  Pulse 
7.  Satisfied  ->  INCLUDE  Current  Detected  Pulse 
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if  TmpWdth  >  PlsMin  7»  Include  Current  Pulse 

°/0  Decrement  While  Loop  Counter  For  EVERY  Detected  Pulse 
WhileCnt=WhileCnt -1 ; 
if  TmpWdth  <  PlsMax 
NumDet  =NumDet  + 1 ; 

TmpDb (NumDet )=SmthZmag(LowDex) -MaxVal ; 

PlsWdthC  NumDet ) = TmpWdth; 

PulseLoc  (NumDet )  =PlsLow  ;  7.  Store  Pis  Location  Index 

PlsDet=TmpZ(PlsLow:PlsHgh) ; 

PlsDur  (NumDet )  =length  (PlsDet  )  ;  7.  Store  Pis  Duration 
PulseVec  =  [PulseVec  ,  PlsDet  ]  ;  °i  Unsorted  Pulse  Vector 
NPlsDet  =  NPlsDet  +  l  ;  7.  Update  Detected  Pulse  Counter 
if  NScr  =  =l  7.  Update  Status  to  Screen  ? 
wait bar ( NP1 sDet/MaxPul .BurstCons , 

[’Burst  Number  ’  ,  num2 st r ( NP1 sDet )  ,  ’  of  ’  ,  .  .  . 
num2 st r ( MaxPul ) ,  ’  Detected.’]) 

end 

end 

end 

TmpZmag  (  PlsLow  :  PlsHgh  )  =MinZ_Db  ;  7»  Remove  Current  Det  Pis 
end  7.  Detection  While  Loop  .  .  .  Detect  Next  Pulse70 
7.  End  Main  While  Loop7« 

if  NScr==l  7.  Update  Waitbar  Screen  Status? 
close  (BurstCons) 
display ([ ’  ’ ] ) 

display([’  A  Total  of  ’, num2  st r ( NumDet )  ,  .  .  . 

’  Pulses  Satisfied  Pulse  Width  Constraints.’]) 
display ([ ’  ’ ] ) 

end 

7.  Put  Detected  Pulses  in  Matrix  Form  with  ONE  pulse  per  row. 
if  NumDet  >  0 

TmpVec  =  PulseVec;  7«  Unsorted  Pulse  Vector 
MaxWidth  =  max(PlsDur); 

PulseMatrix=zeros (NumDet , MaxWidth) ; 
for  k= 1 : NumDet 
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PulseMatrix(k , 1 : PlsDur (k) ) =TmpVec (1 : PlsDur (k) ) ; 

165  TmpVec ( 1 : PlsDur (k) )=[] ; 

end 

7.  Reorder  Pulses  to  Original  Collection  Time  Order 
PlsMat  =  []  ; 

[SortVal , SortLoc] =sort (PulseLoc) ; 

170  PlsMat =zeros (NumDet , MaxWidth) ; 

TmpWdth=PlsWdth ; 
for  k= 1 : NumDet 

PlsMat (k,  :)=PulseMatrix(SortLoc(k)  ,  :)  ; 

PI sWdth (k ) =TmpWdth ( Sort Loc (k) ) ;  "/,  Reorder  Pulse  Widths 

175  PlsDb (k) =TmpDb ( SortLoc (k) )  ;  %  Reorder  Det  Point  Db 

end 

else 

if  NScr==l  %  Update  Waitbar  Screen  Status? 
di splay ( [ 1  ’ ] ) 

180  display([’No  Detected  Pulses  Satisfy  Pulse  Width  ... 

Constraint ’ ] ) 

end 

end 

%Begin  Plotting  Code 

if  NPlot==l  /  Satisfied  ->  Produce  Plots 
185  figure  (1)  %  Magnitude  of  Input  Signal  Plot 

subplot (3,1,1) 
plot ( abs (Z) ) 
grid 

axis  tight 

190  title (’ Magnitude  of  Input  Signal  Z’) 

ylabel (  ’  I Z  I  ’ ) 
xlabel (’ Sample  Number1) 

7. 

subplot (3,1,2) 

195  plot ( SmthZmag) 

grid 

axis  tight 
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title ( [ ’UN -NORMALIZED  Smoothed  |Z|  for  NumSmth  = 
num2str ( NumSmth) ] ) 
ylabel (  1  I Z I  (dB)  ’  ) 
xlabel ( 1  Sample  Number’) 

7. 

subplot (3,1,3) 

plot ( SmthZmag  -  MaxZ_Db) 

grid 

axis  tight 

title ([’NORMALIZED  Magnitude  of  Z  for  NumSmth  =  ’  ,  .  .  . 

num2  st r ( NumSmth) ] ) 
ylabel ( ’ I Z I  (dB) ’ ) 
xlabel (’ Sample  Number’) 

if  NumDet  >  0  7.  Only  Generate  Plots  If  Pulses  Are  Detected 
7,  Create  Sorted  ’VECTOR’  of  Final  Pulses  for  Plotting 
SortVec  =  reshape(PlsMat .  ’  ,1 , NumDet  *M ax Width)  ; 
figure  (2)  7»  Magnitude  of  Input  Signal  Plot 
subplot (3,1,1) 
plot(abs(Z)) 
axis  tight 
grid 

title (’ Magnitude  of  Input  Signal  Z’) 
xlabel (’ Sample  Number’) 
ylabel ( ’  I Z I’  ) 

7. 

subplot (3,1,2) 
plot (abs (PulseVec) ) 
axis  tight 
grid 

title([’ABS  [Unsorted  Pulses]:  ’, num2  st r ( NumDet )  ,  .  .  . 

’  Pulses  Detected’]) 
xlabel (’ Sample  Number’) 
ylabel ( ’  I Z I’  ) 

1 

subplot (3,1,3) 
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plot(abs(SortVec)) 

axis  tight 

grid 

title ([’ABS  [Sorted  Pulses]:  ’  , num2  st r ( NumDet )  ,  .  .  . 

’  Pulses  Detected’]) 
xlabel (’ Sample  Number’) 
ylabel ( ’  I Z |  ’  ) 

7. 

figure  (3)  7.  Relative  Pulse  Amplitude  Plot 
subplot (3,1,1) 
plot (PlsDb  ,  ’  *  ’  ) 
if  Thresh  <  0 

title ([’Rel  Pulse  Amp  at  AddSamp  +  1  = 

num2 st r ( AddSamp  + 1 )  ,  ’  for  Input  Threshold  =  ’  ,  . 

num2 st r ( Thr e sh )  ,  ’  dB  ’  ] ) 

else 

title(’Rel  Pulse  Amp  at  Threshold  Pt :  No  Input 
Threshold ’ ) 

end 

axis  tight 

set(gca,’XLim’,[.98  1. 01*NumDet ] ) 
xlabel (’ Pulse  Number’) 
ylabel ( ’ dB ’ ) 
grid 

7. 

subplot (3,1,2) 
hold 

for  k= 1 : NumDet 

plot(abs(PlsMat(k, :))) ; 

end 

grid 

title ([’Overlay  of  ABS  [PlsMat]: 

num2 str ( NumDet ) , ’  Pulses ,  ’ ,  ... 

’NumSmth  =  ’ , num2 st r ( NumSmth) ] ) 
xlabel (’ Sample  Number’) 
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ylabel ( ’ ABS ’ ) 
axis  tight 

270  subplot (3 , 1 , 3) 

plot(mean(abs(PlsMat)))  ; 
grid 

axis  tight 

title ( [ ’ Column -Wise  Mean  of  ABS  [PlsMat]: 

275  num2str ( NumDet ) , ’  Pulses’]) 

xlabel (’ Sample  Number’) 
ylabel ( ’ Mean ’ ) 

end 

end 

280  /  End  Pulse  Detect  Function 

A.  2  Preamble  Location 

Listing  A. 2:  Code/Locate/LocatePreamble.m 

1  function  [Index , Preamble]  =  LocatePreamble ( Signal , LocMeth , VtThresh ..  . 
, Threshold , SNRdb , Statel , StateQ , IndFlag , Index , F_BW) 
dir  =  ’ F : \ Chamber\ ’ ; 

°/0  Hard  code  to  speed  up  execution 
FilterFreqsHardCode ; 

5  [B,S]  =  size (Signal ) ; 

B. range  =  1 : B ; 

P  =  380; 

Preamble  =  zeros (B ,P)  ; 
if  IndFlag  ~=  1 
10  Index  =  zeros  (B,l); 

end 

Buf  f  er  =  500  ; 

C  =  S  +  Buf  f  er ; 
trans.truth  =  Buffer+1; 

15  RandDat aSt art  =  trans_truth  +  475;  °/0  After  Symbol  region 
wind  =  20 ; 
s  =  2; 
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win_start  =  1  :  s  :  s*floor((  2~11  -wind)/s); 

W  =  length (win_start) ; 

f _RandDat aSt art  =  f ind ( win_ start < RandDat aSt art ,  1,  ’last1  ); 

[Faf,  Fsf]  =  FSfarras; 

[af ,  sf]  =  dualfiltl; 

7»  Set  Chebyshev  type  I  filter  parameters: 

N  =  6;  /  Order  of  filter  --  change  to  6 

Rp  =  0.01;  •/.  Passband  ripple  in  dB 
Tsig  =  XDelta*S; 

Nsig  =  round (Tsig/XDelta) ; 

7.  Set  up  more  filter  parameters 
Fs2  =  l/(2*XDelta) ; 

Wn  =  F_BW/Fs2 ; 

[num  den]  =  chebyl  (N  ,  Rp  ,  Wn)  ;  7«  Filter  coefficients 

for  b  =  B.range  7.  burst 
Data  =  Signal (b , : ) ; 

Data  =  [zeros  (  1  ,  1000)  , Data , zeros  (  1  , 1000) ]  ;  7,  Zero  pad  data 
prior  to  filtering 

Data  =  f  iltf  ilt  (num  ,  den  ,  Data)  ;  7»  Filter  data 
Data  =  Data  ( 1001  :  1000  +  Nsig)  ;  7«  Un-Zero  pad  data  after  ... 
f ilter ing 

Data  =  Data  -  mean(Data); 

Spow  =  sum ( abs ( Dat a) . ~ 2 ) /S ; 

Data  =  Dat a . / ( sqrt ( Spow /2 ) * ( 1+ j ) ) ; 

Spow  =  sum ( abs ( Dat a) . ~ 2 ) /S ; 
nz_rz_I  =  randn(l,C); 
nz_rz_Q  =  randn(l,C); 
nz  =  (nz_rz_I+j *nz_rz_Q) ; 

7.  Filter  Noise 
Nlen  =  length(nz); 

Noise  =  [zeros  (  1  , 1000)  , nz , zeros ( 1  ,  1000) ]  ;  °/„  Zero  pad  noise 
prior  to  filtering 

Noise  =  f  iltf  ilt  (num  ,  den  ,  Noise  )  ;  7»  Filter  noise 
Noise  =  Noise  ( 1001  :  1000  +  Nlen)  ;  7»  Un-Zero  noise  data  after 
f i It  er ing 
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Noise  =  Noi se -mean ( Noi se ) ; 

NP  =  sum ( abs ( No i se )  . ~ 2 ) /C  ; 

Npow  =  ( Spow) / ( 1 0 ~ ( SNRdb / 10) ) ; 
noise  =  sqrt (Npow/NP)  *  Noise; 
sig  =  [zeros (1 , C-S) .Data]  +  noise; 
sig  =  sig ( 1 :  2  ~ 1 1  )  ; 

if  IndFlag  ~=  1 

'/.Feature  Extraction 
windowed  =  zeros (W , wind+1 ) ; 

if  str cmp ( Lo cMeth ,  ’ DenVt  ’  )  I  strcmp ( LocMeth ,  ’  Vt  1  ) 
if  strcmp ( LocMeth ,  1 DenVt  1  ) 
max_level  =  4; 
sig  =  sig-mean ( sig) ; 

y  =  dualtree ( sig , max_level , Faf , af ) ; 
for  p  =  1 : max.level 

aa  =  abs  (y{pMl>)  ; 
bb  =  abs  (y{pH2>)  ; 
cc  =  sqrt((aa).~2  +  (bb).~2); 

Y{p}=y{p>  ;  7. 

7.  Zeroes  coeffs  that  don’t  represent  enough  of... 
value 

[m2,n2]  =  f ind ( abs ( cc ) <Thr e shold (p ) ) ; 

Y{p}{l}(n2)  =  0; 

Y{p}{2}(n2)  =  0; 

end 

Y{p  +  lHl}  =  y{p  +  lHl}; 

Y{p+1}{2}  =  y{p+l}{2}; 

Sig  =  ( idualtree (Y , max.level  , Fsf  , sf ) )  ; 

elseif  str cmp ( LocMeth , ’ Vt 1 ) 
sig  =  sig-mean ( sig) ; 

Sig  =  sig; 

end 

for  w  =  1  : W 

windowed (w , 1 : wind  + 1 )  =  Sig (win.start (w)  : win.start ( . .  . 
w)+wind);  7«  windowing  the  total  Signal 
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=  (Sig(win_start(w) :win_start(w)+wind-l)) ; 
x  =  x-mean (x) ; 

V(w)  =  var(abs(x)); 

85  end 

VT  =  abs (V (1 : end-1) -V (2 : end) ) ; 

[f_index_vt]=DetThresh(VT (1 : f_RandDataStart)  ,  V  (1 :  . .  . 

f _RandDat aSt art ) ,VtThresh) ; 

Index(b)  =  win. s t ar t ( f _ index. vt )  +  wind  -  s; 
elseif  strcmp (LocMeth ,  ’ Fractal  1 ) 

90  sig  =  sig -mean ( s ig) ; 

Sig  =  sig; 
for  w  =  1  :  W 

windowed (w , 1 : wind+ 1 )  =  S ig ( win.st art (w ) : win. st art ( 
w)+wind);  "/«  windowing  the  total  Signal 

end 

95  k  =  1:  wind/2;  °/0  repeat  k  from  1  to  kmax 

fractal.mag  =  Cal cFr act al s ( abs ( windowed) , k) ’ ; 

'/.  Detect  Transient 

[f .index.mag.f rac_pdf , pr ob.mag.f r ac] =Bs cd ( f r act al.mag . 

,3)  ; 

Index(b)  =  win.start (f_index_mag_frac.pdf)  +  wind/2+s; 
100  elseif  strcmp ( LocMeth  ,  1 Perf ’ ) 

Index(b)  =  trans.truth; 
elseif  strcmp ( LocMeth ,  1 Perf Jitter  ’  ) 

Jitter  =  6; 

x  =  round( Jitter  +  ( - Ji tt er - J itt er )  *  rand(l) ) ; 

105  Index(b)  =  trans.truth  +  x; 

end 

end 

Preamble (b ,: )  =  s ig ( Index (b ): Index (b ) +P - 1 ) ; 

end 


Listing  A. 3;  Code/Locate/FilterFreqsHardCode.m 

1  FreqValidMax  =  5 . 189169766750000e+009 ; 

FreqValidMin  =  5  .  170615079250000e  +  009 ; 


107 


XDelta 


4 . 2 105263 1578 9474e- 008 ; 


Fs  =  23750000; 


Listing  A. 4:  Code/Locate/CalcFractals.m 


1  function  [d]  =  CalcFractals (windowed , k) 

'/.Calculates  the  fractal  dimension,  d 
index  =  1 ; 

[M,N]  =  size (windowed) ; 

5  L  =  zer  os  ( length  ( k)  ,  length  (k)  ,  M)  ;  "/0  initialize  for  sum  over  i 

for  a  =  1:  length  (k)  7,  k  new  time  series 

m  =  l:k(a);  %  m  ranges  from  1  to  k 

for  b  =  l:length(m) 

if  (m  (b) +k  (a)  )  <=  (N) 

10  L(b,a,:)  =  sum ( abs ( windowed (:, m (b ) +k ( a)  : k ( a)  : end) . 

windowed/ :  ,m(b)  :k(a)  : end -k ( a) ) )  ,2)  ; 

L (b , a , : )  =  (L (b , a , : ) * (N-l ) / ( f loor ( (N-m (b) ) /k ( a) ) *  k(a) 
))/  k (a) ; 

end 

if  isnan (L (b , a  ,  :  ) ) 
temp  =  0; 

15  end 

end 

end 

L  =  sum(L,l);  7«  average  over  m 
k  =  repmat  (k  ,  [1  ,  1  ,  M]  )  ;  70  repmat  k  to  polyfit 
20  L=squeeze (L) ; 
k=squeeze (k) ; 
for  i  =  1  :  M 

p  =  polyfit(log(k(:,i)),log(L(:,i)),l);  %  least  square  line  fit 
for  log-log 
d ( i )  =  -p  (  1 )  ; 

25  end 


Listing  A. 5:  Code/Locate/Bscd.m 

1  function  [  index , prob_dens ity] =Bs cd ( f r act als , w) 
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N  =  length ( f r act al s ) ; 
prob_density  =  zeros (1 ,N) ; 
m  =  ceil (w/2) ; 

5  N  =  w ; 

for  a  =  1 : N-w 

7,  Piecewise  fractals  for  computer  precision  sake 
d  =  fractals (a : a+w) ; 

p  =  l/(  sqrt  (m*  (N-m)  )  )  /  (  (  sum  (d  .  ~2)  -  (  1/m)  *  sum  (d  (  1  :  m)  )  "2  -  (  1/N- .  .  . 

m) *  sum (d(m  +  l:N))~2)”((N-2)/2))  ; 

10  [prob_density (a+m) ]  =  p; 

end 

[nothing,  index]  =  max (prob_density) ; 

Listing  A. 6:  Code/Locate/Det Thresh. m 

1  function  [ index ] =DetThre sh ( tr a j e ctory , support , threshold) 

N=length ( traj ectory) ; 
w_index  =  0; 

W  =  N  ; 

5  w_index  =  0; 
n=4 . 5 ; 
m=200 ; 

trigger  =  mean ( t raj e ct ory ( 1 : m) ) +n* std ( traj ect or y ( 1 : m) ) ; 
ensure  =  max ( support ( 1 : m) ) ; 

10  trigger  =  threshold ( 1 ) ; 
ensure  =  threshold (2) ; 
for  i  =  1  +  5 :  W 

detect  =  t r a j e ct ory  (  i -5 )  ; 
verify  =  mean ( support ( i -4 : i )) ; 

15  if  detect  >  trigger  &&  verify  >  ensure 

w_index=i -4 ; 
break 

end 

end 

20  if  w_index  ==  0 
w_ index  =  W ; 
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end 


index  =  w_index; 

A. 3  Feature  Extraction 

Listing  A. 7:  Code/Extract /ExtractFeatures.m 

function  Features  =  ExtractFeatures ( Signal , FtrMeth) 

[Faf,  Fsf]  =  FSfarras; 

[af ,  sf]  =  dualfiltl; 

[B,S]  =  size (Signal ) ; 

B_r ange  =  1 : B ; 
switch  FtrMeth 
case  ’ Td ’ 

Features  =  zeros (B  ,  3 , 9)  ; 
case  { ’ Wd ’ } 

Features  =  zeros (B  ,  3 , 5 , 9)  ; 

end 

Region  =  { ’ pre 1 1 , ’ pre2 ’ , ’ all ’ } ; 
trans_length . pre 1  =  190; 
trans.length . pre2  =  190; 
trans_length . all  =  380; 
for  b  =  B.range  %  burst 
sig  =  Signal (b  ,  :  )  ; 
for  x  =  1 : length ( Region ) 

if  x  ==  2  '/,pre2  -  must  bypass  all  of  prel 
index.start  =  l+trans_length . prel ; 

else 

index.start  =  1; 

end 

index.end  =  index_start+trans_length . (Region{x})  -  1; 
Sig  =  sig ( index.start : index.end) ; 

Sig  =  Sig  -  mean ( Sig) ; 
switch  FtrMeth 
case  1  Td  ’ 

Fe atur e s (b , x , : )  =  InstFeatures ( Sig) ; 
case  { ’ Wd ’ } 
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end 


end 

end 


zp  =  2 . ~ ce il ( log2 ( length ( Sig) ) )  -  length(Sig); 

switch  FtrMeth 
case  ’ Wd ’ 

SIG  =  real (Sig) ; 

end 

Coeffs  =  dualtree ( [SIG ,  zeros (1 ,zp) ] ,max_level ,Faf 
,af); 

for  p  =  1 : max_level+l 

if  p  ==  max_level+l 

IndTemp  =  ce il ( length ( SIG) /( 2 ~ max.level )) ; 

else 

IndTemp  =  ce il ( length ( SIG) /( 2 ~ p )) ; 

end 

switch  FtrMeth 

case  { ’ Wd ’ , ’ WdC ’ > 

aa  =  Coef f s Ip }{ 1 }( 1 : IndTemp) ; 
bb  =  Coef f s Ip }{2} ( 1 : IndTemp) ; 

Feat ur e s (b , x , p , : )  =  InstFeatures (aa  + 

j  *bb )  ; 

end 

end 


Listing  A. 8:  Code/Extract/InstFeatures.m 

1  function  [ Ins t .Features ]= Inst Feat ures ( s ignal ) 
i  =  [ 1 : length ( s ignal ) ] 

Fs  =  23.T5E6; 

Tsamp  =  1/Fs ; 

51  =  real ( signal ) ; 

Q  =  imag (signal ) ; 

Unwrap_Phase  =  unwr ap ( at an2 ( Q , I ) ’ ) ; 

Inst_Freq  =  gr adi ent ( Unwr ap_Phase , Tsamp ) / (2* pi ) ’ ; 
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Inst_Amp  =  abs ( signal )’ ; 

10  iu_f  =  mean ( Inst_Freq) ; 

Inst_Phase  =  Unwrap_Phase  -  2*pi * i *mu_f /Fs ; 

Inst_Freq  =  Inst_Freq  -  mu_f ; 

Inst_Phase  =  Inst_Phase  -  mean ( Inst_Phase) ; 

Inst_Amp  =  Inst_Amp  -  mean ( Inst_Amp) ; 

15  “/«  Calculate  variance  of  sub-segment 
Amp.var  =  var ( Inst_Amp) ; 

Pha_var  =  var ( Inst.Phase) ; 

Fre.var  =  var ( Inst_Freq) ; 

7,  Calculate  skewness  of  sub  -  segment 
20  Amp_skew  =  skewness ( Inst_Amp) ; 

Pha_skew  =  skewness ( Inst.Phase) ; 

Fre_skew  =  skewness ( Inst_Freq) ; 

7,  Calculate  kurtosis  of  sub  -  segment 
Amp_kurtosis  =  kurtosis ( Inst_Amp) ; 

25  Pha.kurtosis  =  kurtosis ( Inst_Phase) ; 

Fre.kurtosis  =  kurtosis ( Inst_Freq) ; 

Inst.Features  =  [Amp.var  ,  Pha_var  ,  Fre.var  ,  Amp_skew  ,  Pha_skew , 
Fre_skew ,  Amp_kurtosis ,  Pha.kurtosis ,  Fre.kurtosis] ; 

A. 4  Device  Classification 

Listing  A. 9:  Code/Classify/ClassifyDevices.m 

1  function  ClassAcc  =  Classif yDevices (Dvcl , Dvc2 , Dvc3) 

[B , dim]  =  size(Dvcl); 
k_f  old  =  5; 
if  k_fold  ==  1 
5  n.class  =  1:B; 

n_tr ain  =  1 : B ; 

else 

n.class  =  1 : ceil ( l/k_f old*B ) ; 
n.train  =  n.class ( end) +1 : B ; 

10  end 

N_class  =  length (n.class) ; 

7,  Initialize  Confusion  Matrix  Variables  within  the  SNR  loop 
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AA  =  0; 

AB  =  AA; 

15  AC  =  AA; 
Aerr  =  AA  ; 


BB  = 

AA; 

BA  = 

AA; 

BC  = 

AA; 

20  Berr 

=  AA; 

CC  = 

AA; 

CA  = 

AA; 

CB  = 

AA; 

Cerr 

=  AA; 

classldata  = 
class2data  = 
class3data  = 
7.  %  Relevant 


reshape (Dvcl , [B,dim] ) 
reshape(Dvc2 , [B , dim] ) 
reshape(Dvc3 , [B , dim] ) 
Features 


30  7.  f  27  =[102, 10, 43, 57,91, 135, 97, 110,114, 24, 33, 7, 65, 105,12, 99, 132, 
7.31  ,92 ,30 , 14 ,4 , 130 , 18 ,39 , 17 , 108]  ; 

7.  classldata  =  classldata(:,f27); 

7.  class2data  =  class2data(:,f27); 

7.  class3data  =  class3data(:,f27); 

35  classltrain  =  classldata; 
class2train  =  class2data; 
class3train  =  class3data; 

7.  Scramble  data 
scramble  =  randperm(B); 

40  classldata  =  classldata] scramble ,:) ; 
class2data  =  class2data( scramble ,:) ; 
class3data  =  class3data( scramble ,:) ; 
classltrain  =  clas s It rain ( s cr amble ,:) ; 
class2train  =  class2train ( scramble ,:) ; 

45  class3train  =  class3train ( scramble ,:) ; 
ave_k  =  zeros  ( 1  ,  k_f  old  )  ; 
for  k  =  1 : k_f old 
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i  =  n.class; 

classldata_k  =  classldata(i , : ) ; 
class2data_k  =  class2data(i , : ) ; 
class3data_k  =  class3data(i , : ) ; 
i  =  n_train; 

class ltrain_k  =  class ltrain ( i ; 
class2train_k  =  clas s2t r ain ( i , : ) ; 
class3train_k  =  clas s3t rain ( i ,:) ; 

classldata  =  c ir c shif t ( class ldat a ,  [N_class  ,  0] )  ; 
class2data  =  c ir c shif t ( clas s2dat a , [N_clas s , 0] ) ; 
class3data  =  c ir c shif t ( clas s3dat a  ,  [N_clas s  ,  0] )  ; 
classltrain  =  circshif t ( classltrain ,  [N_class  , 0] )  ; 
class2train  =  circshif t ( class2train , [N_class , 0] ) ; 
class3train  =  circshif t ( class3train  ,  [N_class  , 0] )  ; 

7.  Training  the  MDA  Projection  Matrix  and  ML  paramters 
[XI ,X2 ,xl ,x2 ,R1 ,R2 ,R3 , data.mean , data_std ,W, Fishclassl  , . .  . 
Fishclass2  , Fishclass3] =Train_class2 ( classltrain_k , . .  . 
class2train_k , class3train_k , classldata_k , class2data_k  , . .  . 
class3data_k) ; 

•/.  Clas  s  i  f  y  ing 

[Tot_err , N_tot , A_in_A , A_in_B , A_in_C , A_err , B_in_B , B_in_A , B_in_C 
, B_err , C_in_C , C_in_A , C_in_B , C_err ,  Fishclassl,  Fishclass2, 


F ishclass3] . . . 

=  Classify_class(classldata_k  , class2data_k  , class3data_k  ,X1 , 
X2 ,xl ,x2 ,R1 ,R2 ,R3 , data.mean , data_std ,  W) ; 
ave_k(k)  =  (N_tot -Tot.err) . /N_tot  ; 

7,  7.  Confusion  Matrix 


7, 

7. 

7. 

7. 

7. 

7. 

7. 

7. 


AA  =  AA  + A_in_A  ; 

AB  =  AB+A_in_B ; 

AC  =  AC+A_in_C ; 
Aerr  =  Aerr+A_err; 
BB  =  BB+B_in_B ; 

BA  =  BA+B_in_A ; 

BC  =  BC+B_in_C ; 
Berr  =  Berr+B_err; 
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7.  CC  =  CC+  C_  in_C 

7,  CA  =  CA+  C_  in_  A 

7.  CB  =  CB  +  C_  in_B 


80  7. 


Cerr  =  Cerr+C_err; 


end 

ClassAcc  =  ave_k; 


Listing  A.  10:  Code/Classify/Train. m 


1  function  [XI  X2  xl  x2  R1  R2  R3  data_mean  data_std,  W,  Fishclassl , 
Fishclass2  ,  Fishclass3]=Train_class2(classldata , class2data , . .  . 
class3data ,gridl ,grid2 ,grid3) 


Limit 

=  3; 

Res 

=  .01; 

ResPts  = 

400  ; 

Dis 

=  2; 

N_rec 

=  length ( c las s ldat a ( :  ,1 

N_tot 

=  N_rec*3; 

N_recg 

=  length ( gridl (:,  1 ) )  ; 

N_totg 

=  N_recg  *3 ; 

10  1111111111 

1  PART  I  7. 

7  7  7  7  7  7  7  7  7  7 

/o  /o  /0  / 0  /o  /o  /o  /o  /o  /o 


7  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  / 0  /o 


7.  Globally  normalize 
15  training.data  = 

data.mean 
data_std 

training_data_norm  = 
classldata 
20  class2data 
class3data 


emission  records 

[classldata; class2data; class3data] ; 
ones (N_tot , 1 ) *mean ( t r aining_dat a) ; 
ones (N_tot ,l)*std(training_data,l) ; 
(training_data-data_mean) . / data_std; 
training_data_norm(N_rec*0+l : N_rec*l 
training_data_norm(N_rec*l+l : N_rec*2 
training_data_norm(N_rec*2+l : N_rec*3 


11111111111 


) 

) 

) 


1  PART  II  1 


7  7  7  7  7  7  7  7  7  7  7 

/o  /o  /o  /o  /o  /o  / 0  /o  /o  /o  /o 
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25  °/o  Compute  the  within  class  covariance  matrixes  S_k  ,  where  k  = 
1,2,3 

51  =  cov ( class ldata) ; 

52  =  cov ( class2data) ; 

53  =  cov ( class3data) ; 

7,  Compute  the  SW  matrix  by  summing  the  S_k  matrixes 
30  SW  =  SI  +  S2  +  S3  ; 


7,  Compute  the  within  class  and  total  means 
ml  =  mean ( class ldata) 
m2  =  mean ( class2data) 
m3  =  mean ( class3data) 

35  mtot  =  ( 1 / N_t ot ) * ( N_r ec *ml  +  N_rec*m2  +  N_rec*m3); 

7.  Compute  the  SB  matrix 

SB  =  N_rec*(ml  -  mtot) * ( (ml  -  mtot) ’)  +. . . 

N_rec*(m2  -  mtot)*((m2  -  mtot)  ’)  +.  .  . 

N_rec*(m3  -  mtot)*((m3  -  mtot)1); 

40  7«  Solve  for  x 
x  =  SW\ SB ; 

7.  Find  Fisher  plane  matrix  W  (eig  vectors  corresponding  to  the  two 
largest  eig  values) 

[V , D]  =  eig(x) ; 

lamda  =  sum ( D ) ; 

45  maxi  =  find(lamda  ==  max(lamda)); 

lamda(maxl)  =  NaN; 

max2  =  find(lamda  ==  max(lamda)); 

W  =  [V(:  ,maxl)  ’;V(:  ,max2)  ’]  ; 

7,  Normalize  each  emission  record  using  the  parameters  calculated  . 
during  training 

50  training_datag  =  [gridl ; grid2 ; grid3] ; 

data.meang  =  ones (N_totg , 1) *data_mean (1 , : ) ; 

data_stdg  =  ones (N_totg , 1) *data_std (1 , : ) ; 

training_data_normg  =  ( t r aining_dat ag - dat a.meang) . / dat a_ st dg ; 


classldatag 


=  training_data_normg(N_recg*0+l : N_recg*l  , 
=  training_data_normg(N_recg*l+l : N_recg*2 , 
=  training_data_normg(N_recg*2+l : N_recg*3 , 


)  ; 
)  ; 
)  ; 


55  class2datag 
class3datag 
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Fishclasslg  =  W  *  class ldatag  ’  ; 

Fishclass2g  =  W  *  class2datag ’ ; 

Fishclass3g  =  W  *  class3datag  ’  ; 

7.7.7.  7. 7. 7. 7. 7. 7. 7. 7, 

7.  PART  II  7, 

7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 

7«  Project  each  class  using  the  Fisher  plane  calculated  during  ... 
training 

Fishclassl  =  W  *  classldata ' ; 

Fishclass2  =  W  *  class2data ' ; 

Fishclass3  =  W  *  class3data  '  ; 

LoLimitl  =  min([Fishclasslg(l,:) ,Fishclass2g(l,:) ,Fishclass3g(l,: 

, Fishclassl (1, :) ,Fishclass2(l, :) ,Fishclass3(l , :)]) ; 

LoLimit2  =  min([Fishclasslg(2,  :)  ,Fishclass2g(2,  :)  ,Fishclass3g(2,  : 

,Fishclassl(2, :) ,Fishclass2(2, :) , Fishclass3(2, :)]) ; 

HiLimitl  =  max([Fishclasslg(l,:) ,Fishclass2g(l,:) ,Fishclass3g(l,: 

, Fishclassl (1, :) ,Fishclass2(l, :) ,Fishclass3(l , :)]) ; 

HiLimit2  =  max([Fishclasslg(2, :) ,Fishclass2g(2, :) ,Fishclass3g(2, : 

,Fishclassl(2, :) ,Fishclass2(2, :) ,Fishclass3(2, :)]) ; 

ResPts  =  400 ; 

xl  =  1 inspac e ( LoLimi t 1 , HiLimit 1 , ResPt s ) ; 

x2  =  linspace ( LoLimit2 , HiLimit2 , ResPts ) ; 

77.7.7.7.7.7.7.7.7.7.7. 

7.  PART  III  7. 

7  7  7  7  7  7  7  7  7  7  7  7 

/o  /o  / 0  /o  /o  /o  /o  /o  /o  /o  /o  / 0 

7.  Find  the  mean  vector  for  each  class 
meanl  =  mean (Fishclassl 1 ) ; 
mean2  =  mean (Fishclass2 1 ) ; 
mean3  =  mean (Fishclass3  1  )  ; 

7.  Find  the  covariance  matrix  for  each  class 
K1  =  cov ( Fi shclas s 1  ’)  ; 

K2  =  cov ( Fi shclas s2  ’)  ; 

K3  =  cov ( Fi shclas s3  ’)  ; 

7.  Find  the  inverted  covariance  matrix  for  each  class 
Q1  =  inv (K1 )  ; 
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Q2  =  inv ( K2 )  ; 

Q3  =  inv ( K3 ) ; 

[XI  X2]  =  meshgr id (xl , x2 )  ; 

90  Kldet  =  (  1/ (2* pi  * sqrt ( det (K1 )  )  )  )  ; 

K2det  =  ( 1 / ( 2 *  pi  * sqrt ( det ( K2 ) ) ) )  ; 

K3det  =  ( 1 / ( 2* pi  * sqrt ( det ( K3 ) ) ) )  ; 
pxl  =  zeros ( length (x2 ), length (xl ))  ; 
px2  =  pxl ; 

95  px3  =  pxl ; 

for  i  =  l:length(xl) 

for  k  =  l:length(x2) 

pxl(k,i)  =  [ ( xl ( i ) -meanl ( 1 ) ) , ( x2 (k) -meanl ( 2) ) ] * Q 1 

* [(xl (i) -meanl (1)) ; (x2(k)-meanl (2))] ; 
px2(k,i)  =  [ ( xl ( i ) -mean2 ( 1 ) ) , ( x2 (k) -mean2 ( 2) ) ] * Q2 

*[(xl(i)-mean2(l)) ; (x2(k)-mean2(2))] ; 

100  px3(k,i)  =  [  ( xl ( i ) -mean3 ( 1 ) )  , ( x2 (k) -mean3 ( 2) ) ] * Q3 

*[(xl(i)-mean3(l)) ; (x2(k)-mean3(2))] ; 

end 


end 

pxl  =  Kldet 

*  exp (  - 

•  5*(pxl) )  ; 

px2  =  K2det 

*  exp (  - 

5* (px2 ) )  ; 

105  px3  =  K3det 

*  exp  (  - 

5* (px3) )  ; 

7,  Initialize  the  regions 
R1  =  zeros ( size (pxl ))  ; 

R2  =  Rl; 

R3  =  Rl; 

110  i  Define  the  Bayesian  decision  regions 
Rl(pxl>=px2  &  pxl>px3)  =  1; 

R2(px2>pxl  &  px2>=px3)  =  1; 

R3(px3>=pxl  &  px3>px2)  =  1; 

[DD,LL]  =  bwdist (R1+R2+R3) ; 

115  Rl (f ind(Rl (LL) ) )  =  1 ; 

R2(find(R2(LL)))  =  1; 

R3(find(R3(LL) ) )  =  1; 
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7«  figure 

120  i  hold  on 

7.  colormap( [1  0  0 ; 0  1  0 ; 0  .5  1]) 


7.  plot  (-10, 

O 

T — 1 

1 

’s’, 

, ’ MarkerFaceColor ’ 

,  'r  ’  , 

, ’ MarkerEdgeColor ’ 

,’r 

’) 

/  plot (-10, 

1 

h- »■ 

O 

’  s  3  , 

, ’MarkerFaceColor’ 

,’g’: 

, ’MarkerEdgeColor’ 

,’g 

’) 

7.  plot  (-10, 

o 

T — 1 

1 

}  s  3  . 

, ’MarkerFaceColor’ 

,’b’: 

, ’MarkerEdgeColor’ 

,  ’b 

’) 

125  7.  plot  (-10  , 

o 

T — 1 

1 

J  > 

,  ’MarkerFaceColor ’ 

,’k’: 

, ’MarkerEdgeColor’ 

,  ’k: 

’) 

7«  piot(-io, 

1 

1-^ 

O 

^  0  J  . 

, ’MarkerFaceColor’ 

,’k’: 

, ’MarkerEdgeColor’ 

,  ’k: 

’) 

7.  Floor_level  =  -(max([max(max(pxl))  max(max(px2))  max(max(px3))]) 

)  ; 

7«  7» Combine  and  plot  3D  Gaussians  with  mean  and  covariance  equal  to 
the  class  mean  and  covariance 
70  px_sum  =  pxl.*Rl  +  px2.*R2  +  px3.*R3; 

130  /  Gousians  =  surf (XI  , X2 , px_sum ,’ LineStyle ’,  ’ none  1  )  ; 

/  for  i  =  l:(min(size(Xl))-l)/24:min(size(Xl))-l 
7«  i  =  round  (  i  )  ; 

7,  plot3(Xl(:  ,i)  ,X2(:  ,i)  ,  px_sum  ( :  ,  i )  ,  ’  k  ’ )  ; 

7.  plot3(Xl(i,:),X2(i,:),px_sum(i,:),’k’); 

135  7.  end 

7.  set  (Gousians  ,  ’  Cdatamapping  ’  ,  ’direct  ’  ) 

7.  set  (Gousians  ,  ’Cdata’  , 1*R1+2*R2+3*R3) 

7.  7.Plot  the  projected  points  from  each  class 
7.  scatter3(Fishclassl(l,l:Dis:end)  ,Fishclassl(2,l:Dis:end) 
Floor_level*ones(l, ceil (size (Fishclassl  ,2)/Dis))  ,’r.  ’  ) 

140  7.  scatter3(Fishclass2(l,l:Dis:end)  , Fishclass2(2,l:Dis:end) 
Floor_level*ones(l, ceil (size (Fishclass2  ,2)/Dis))  ,’g.  ’) 

7,  scatter3(Fishclass3(l,l:Dis:end)  , Fishclass3(2,l:Dis:end) 
Floor_level*ones(l, ceil (size (Fishclass3 ,2)/Dis)) ,’b. ’) 

7.  7»Plot  the  means  of  the  projected  points  from  each  class 
7,  scatter3(meanl  (1)  ,meanl  (2)  ,Floor_level  ,  ’ko  ’  ,  ’filled’) 

7.  scatter3(mean2(l)  ,mean2(2)  ,Floor_level  ,  ’ko’  ,  ’filled’) 

145  7.  scatter3  (mean3  (  1 )  ,  mean3  (2)  ,  Floor_level  ,  ’  ko  ’  ,  ’  f  illed  ’  ) 

7,  7. Plot  the  Bayesian  decision  regions  for  the  classes 
l  contour3(Xl,X2,(Rl  +  Floor_level-.5),l,’r’) 

7«  contour3(Xl,X2,(R2  +  Floor_level-.5),l,’g’) 
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7.  contour3(Xl  ,  X2  ,  (R3  +  Floor.level  -  .  5)  ,  1  ,  ’b  ’  ) 

150  i  7, Set  axis  parameters 
7.  xlabel  (  ’  Y1  ’  ) 

7.  ylabel  ( ’ Y2  ’  ) 

7o  legend  (  1  Class  A’, ’Class  B’, ’Class  C’,  ’Class  Point  ’,’  Class  Mean’) 
7.  axis  (  [min(xl)  max(xl)  min(x2)  max(x2)  Floor.level  -Floor.level]) 
155  7.  view  (45, 22. 5) 

7«  lightangle  (45 , 22 . 5) 

7c  light  (’ Style inf  inite ’)  ; 

7c  material  shiny 
7c  lighting  phong 
160  7c  grid  on 

7c  hold  off 


Listing  A. 11:  Code/ Classify /Classify. m 


1  function  [Tot.err , N_tot  , A_in_A , A_in_B , A_in_C , A_err , B_in_B , B_in_A  ,  . 
B_in_C , B_err , C_in_C , C_in_A , C_in_B , C_err ,  Fishclassl,  Fishclass2 
,  F ishclass3] = . . . 

Classify.class (classldata , class2data , class3data ,X1  ,X2  ,xl  ,x2 ,R1 
,R2 ,R3 , dat a .mean , data.std , W) 

Di  s  =  1 ; 

N_rec  =  length ( c las s ldat a ( :  ,1))  ; 

5  N_tot  =  N_rec*3; 

/ 7. 1 7. 7. 7. 7. 7. 7. 7. 

7.  PART  I  7. 


mmrm 

7c  Normalize  each  emission  record  using  the  parameters  calculated 
during  training 

10  training.data  =  [ c 1  as s 1 dat a ; class2dat a ; c las s3dat a]  ; 

data.mean  =  ones (N_tot , 1 ) *data_mean ( 1 , : ) ; 

data.std  =  one s ( N_t ot , 1 ) * dat a_ st d ( 1 , : ) ; 

training_data_norm  =  ( training.dat a - dat a.mean)  ./ dat a_std ; 
classldata  =  training_data_norm (N_rec *0  +  1 : N_rec *1  ,  : ) 

15  class2data  =  t raining.dat a.norm ( N_rec *  1  + 1 : N_re c *2  ,: ) 

=  training_data_norm(N_rec*2  +  l : N_rec*3  ,  : ) 


class3data 
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it/: nr  Lira 

1  PART  II  1 

y  y  y  y  y  y  y  y  y  y  y 

/o  /o  /o  / 0  /o  /o  /o  /o  /o  /o  /o 

7o  Project  each  class  using  the  Fisher  plane  calculated  during 
training 

Fishclassl  =  W  *  classldata ' ; 

Fishclass2  =  W  *  class2data ' ; 

Fishclass3  =  W  *  class3data  '  ; 

1  111111111111 
1  1  PART  III  1 
1  1 1 1 1 1 1 1 1 1 1 1 1 

1  Round  and  rescale  data  for  confusion  matrix  calculations 
scalel  =  size (R1 , 2) -1 ; 
scale2  =  size (R1 , 1) -1 ; 

Xl_shift  =  min (min (XI ) ) ; 

X2_shift  =  min ( min ( X2 ) )  ; 

Xl.scale  =  max(max(Xl))  -  Xl_shift; 


X2_scale  =  max(max(X2) )  -  X2_shift: 


Xl_l 


1  +  r ound ( ( ( Fi shclas s 1 ( 1 , : ) 


scalel 


) ; 


X2  1 


1  +  r ound ((( Fi shclas s 1 (2 ) 


scale2 


) ; 


XI  2 


1  +  r ound ((( Fi shclas s2 ( 1  ,: ) 


Xl_shift)  /  Xl.scale  )  *  ... 


X2_shift)  /  X2_scale  )  * 


Xl_shift)  /  Xl.scale  )  * 


scalel 


) ; 


X2_2  =  1  +  r ound ((( Fi shclas s2 (2 ,: )  -  X2_shift)  /  X2_scale  )  *  ... 


scale2 


) ; 


Xl_3  =  1  +  r ound ((( Fi shclas s3 ( 1 ,: )  -  Xl_shift)  /  Xl.scale  )  *  ... 


scalel 


) ; 


X2_3  =  1  +  r ound ((( Fi shclas s3 (2  ,: )  -  X2_shift)  /  X2_scale  )  * 


scale2 


) ; 


1  Initialize  the  individual  classification  terms 
A  in  A  =  0 ; 


A_ in_B 
A_ in_C 
A_err 
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B_in_A  =  0; 

B_in_B  =  0; 

B_in_C  =  0; 

B_err  =  0; 

C_ in_ A  =  0; 

C_in_B  =  0; 

C_ in_C  =  0; 

C.err  =  0; 

n=  f ind ( X2_ 1  >  0  &  Xl_l  >  0  &  X2_l  <  scale2 
A_in_A  =  sum ( diag (R1 ( X2_ 1 (n) , Xl_ 1 (n) ) ) ) ; 
A_in_B  =  sum ( diag ( R2 ( X2_ 1 (n)  , X 1 _ 1 (n) )))  ; 
A_in_C  =  sum ( diag (R3 ( X2_ 1 (n)  , Xl_  1  (n) )))  ; 
A_err  =  N_rec - length (n) ; 

n=f ind (X2_2  >  0  &  Xl_2  >  0  &  X2_2  <  scale2 
B_in_A  =  sum ( di ag ( R1 ( X2_2 (n)  , X 1 _2 (n) ) ) )  ; 

B_ in_B  =  sum (diag(R2(X2_2(n)  ,Xl_2(n))))  ; 
B_in_C  =  sum ( di ag ( R3 ( X2_2  (n)  , X 1 _2 (n) ) )  )  ; 
B_err  =  N_rec - length (n)  ; 

n=f ind (X2_3  >  0  &  Xl_3  >  0  &  X2_3  <  scale2 
C_in_A  =  sum ( di ag ( R1 ( X2_3 (n) , X 1 _3 (n) ) ) ) ; 
C_in_B  =  sum ( di ag ( R2 ( X2_3 (n) , X 1 _3 (n) ) ) ) ; 
C_in_C  =  sum ( diag ( R3 ( X2_3 (n) , X 1 _3 (n) ))) ; 
C_err  =  N_rec - length (n)  ; 
i  Count  total  misclassifications 
Tot.err  =  A_in_B  +  A_in_C  +  A_err  +  B_in_A 
C_in_A  +  C_in_B  +  C.err ; 

7.  figure 

7.  XI  =  1  +  ((XI  -  Xl_ shift)  /  Xl.scale  ) 
7.  X2  =  1  +  ( ( X2  -  X2_shif  t)  /  X2_scale  ) 
7.  hold  on 

7.  plot  (-10, -10,  ’  x  ’  ,  ’  MarkerFaceColor  ’  ,  ’  r  ’  ,  ’ 
°i  plot  (-10, -10,  ’  +  ’  ,  ’MarkerFaceColor’  ,  ’g’  ,  ’ 
i  plot (-10, -10,  ’ * ’  ,  ’MarkerFaceColor’  ,  ’ b ’  ,  ’ 
i  %  Plot  the  Baysian  decision  regions 
i  contour(Xl ,X2 ,R1 ,1 , ’ r ’ ) 


&  Xl_l  <  scalel ) ; 


&  Xl_2  <  scalel ) ; 


&  Xl_3  <  scalel ) ; 


+  B_in_C  +  B_err  +  . 

*  scalel  ; 

*  scale2  ; 

MarkerEdgeColor ’  ,  ’ r  ’) 
MarkerEdgeColor’  ,  ’g’) 
MarkerEdgeColor’  ,  ’ b  ’  ) 
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7.  contour(Xl  ,  X2  ,  R2 , 1  ,  ’g  ’  ) 

80  7,  contour(Xl  ,X2  ,R3 ,  1  ,  ’  b  ’  ) 

7.  7.  Plot  the  test  points 

7.  scatter(Xl_l(l:Dis:end)  ,X2_l(l:Dis:end)  ,  ’  x  ’  ,  ’  r  ’  ) 

I  scatter(Xl_2(l:Dis:end)  ,X2_2(l:Dis:end)  ,  ’  +  ’  ,  ’  g  ’  ) 

7.  scatter(Xl_3(l:Dis:end)  ,X2_3(l:Dis:end)  ,  ’  *  ’  ,  ’  b  ’  ) 

85  7.  7.  Set  axis  paramiters 
I  xlabel ( ’ Y1 ’ ) 

/  ylabel ( 1 Y2 ’ ) 

7.  legend  (’  Class  A  Point’,1 Class  B  Point’,  ’Class  C  Point’) 

7»  axis([l  scalel  1  scale2]) 

90  7.  s  e  t  (  g c  a  ,  ’  XT  i  c k  ’  ,  (  1  :  (  s  c  a  1  e  1  -  1 )  /  4  :  s  c  a  1  e  1  )  ) 

7.  set(gca,’YTick’  ,  (I:(scale2-l)/4:scale2)) 

7o  set  (gca  ,  ’XTickLabel  ’  ,{min(xl)  ,  ((min(xl)  )  +  (min(xl)+max(xl)  )  /  2  )  .  . 

II,... 

7.  (min(xl)+max(xl))/2,  (((min(xl)+max(xl))/2)+max(xl))/2,  max(xl 

)}) 

7.  set(gca,  ’YTickLabel  ’  ,{min(x2)  ,  ((min(x2))+(min(x2)+max(x2))/2  )  .  . 
/2,  .  .  . 

95  7.  (min  (  x2  ) +max  (  x2  )  ) /2  ,  (  (  (min  (  x2  ) +max  (  x2  )  ) /2) +max  (  x2  )  ) /2  ,  max(x2 

)}) 

7.  axis  square 
7.  grid  on 
7»  hold  off 

A 5  DT-CWT  Transformation 

Listing  A. 12:  Code/DualTree/dualtree.m 

1  function  w  =  dualtree(x,  J,  Faf,  af ) 

7.  Dual-tree  Complex  Discrete  Wavelet  Transform 
7.  USAGE: 

5  7«  w  =  dualtree(x,  J,  Faf,  af ) 

7.  INPUT: 

7.  x  -  N-point  vector 
7o  1)  N  is  divisible  by  2~J 
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7.  2)  N  >=  2  ~  ( J  —  1 )  *length  (af) 

7.  J  -  number  of  stages 
7.  Faf  -  filters  for  the  first  stage 
°i  af  -  filters  for  the  remaining  stages 
7,  OUTPUT  : 

7«  w  -  DWT  coefficients 
7.  w{  j  Hi}  ,  j  =  1..J  -  real  part 

7.  w{jH2},  j  =  1..J  -  imaginary  part 

l  w{J+l}{d}  -  lowpass  coefficients,  d  =  1,2 

7.  EXAMPLE: 

7.  x  =  rand(l,  512); 

7.  J  =  4; 

7.  [Faf,  Fsf]  =  FSfarras; 

7.  [af  ,  sf]  =  dualfiltl; 

l  w  =  dualtree(x,  J,  Faf,  af ) ; 

7.  y  =  idualtree(w,  J,  Fsf,  sf); 

7.  e  r  r  =  x  -  y  ; 

7«  max  (  abs  (  err  )  ) 

7.  WAVELET  SOFTWARE  AT  POLYTECHNIC  UNIVERSITY,  BROOKLYN, 
7.  http  :  /  / taco  .  poly  .  edu/WaveletSoftware/ 

7«  normalization 
x  =  x/ sqrt (2) ; 

7.  Tree  1 

[xl  w{lHl>]  =  afbDT(x,  Faf{l}); 
for  j  =  2 : J 

[xl  w{j  Hi}]  =  af  bDT  (xl  ,  af{l}); 

end 

w{J  +  lHl}  =  xl; 

7.  Tree  2 

[x2  w{l}{2}]  =  af bDT (x ,  Faf{2}); 
for  j  =  2 : J 

[x2  w{j}{2}]  =  afbDT(x2,  af{2}); 

end 

w{  J  +  1H2}  =  x2; 
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Listing  A. 13:  Code/DualTree/afbDT.m 


function  [lo ,  hi]  =  afbDT(x,  af ) 

7o  Analysis  filter  bank 
°/„  USAGE: 


[lo  , 

hi]  = 

af  b  (x  ,  af  ) 

INPUT  : 

x  - 

N-point 

vector ,  where 

1) 

N  is  even 

2) 

N  >=  length ( af ) 

af  - 

analy s 

is  filters 

af  (  : 

,  1)  - 

lowpass  filter 

(even  length) 

af  (  : 

,  2)  - 

highpass  filter 

(even  length) 

OUTPUT  : 

1  lo  -  Low  frequecy  output 
“/.  hi  -  High  frequency  output 
“/„  EXAMPLE: 

7»  [af  ,  sf]  =  farras; 

“/.  x  =  rand  (1 ,64)  ; 

“/o  [lo  ,  hi]  =  afb(x,  af )  ; 

V.  y  =  sfb(lo,  hi,  sf); 

7«  err  =  x  -  y; 

“/o  max  (  abs  (  err  )  ) 

7,  WAVELET  SOFTWARE  AT  POLYTECHNIC  UNIVERSITY,  BROOKLYN, 
i  http : // taco . poly . edu/WaveletSoftware/ 

N  =  length (x)  ; 

L  =  length ( af ) /2  ; 
x  =  cshif t (x , -L)  ; 

I  lowpass  filter 

lo  =  upfirdn(x,  af ( : , 1) ,  1,  2); 

lo ( 1 : L )  =  lo (N/2+ [1 : L] )  +  lo(l:L); 
lo  =  lo  ( 1  :  N/2 )  ; 

%  highpass  filter 

hi  =  upfirdn(x,  af ( : , 2) ,  1,  2); 
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35  hi ( 1 : L)  =  hi (N/2+  [1 : L] )  +  hi(l:L); 
hi  =  hi (1 : N/2)  ; 

Listing  A.  14:  Code/DualTree/idualtree.m 

1  function  y  =  idualtree(w,  J,  Fsf,  sf) 

7.  Inverse  Dual-tree  Complex  DWT 
7.  USAGE: 

5  70  y  =  idualtree(w,  J,  Fsf,  sf) 

7.  INPUT  : 

7»  w  -  DWT  coefficients 

7.  J  -  number  of  stages 

7.  Fsf  -  synthesis  filters  for  the  last  stage 

10  7.  sf  -  synthesis  filters  for  preceeding  stages 
7.  OUTUT  : 

7.  y  -  output  signal 

7,  See  dualtree 

7.  WAVELET  SOFTWARE  AT  POLYTECHNIC  UNIVERSITY,  BROOKLYN, 
15  7.  http://taco.poly.edu/WaveletSoftware/ 

7.  Tree  1 

yl  =  w{J+1}{1} ; 
for  j  =  J : - 1 : 2 

20  yl  =  sf bDT (yl ,  w{j}{l},  sf{l>); 

end 

yl  =  sf  bDT  (yl  ,  w{lHl},  Fsf{l>); 

7.  Tree  2 

y2  =  w{J+l}{2} ; 

25  for  j  =  J : -1 : 2 

y2  =  sf bDT (y2 ,  w{j}{2>,  sf{2>); 

end 

y2  =  sf bDT ( y2 ,  w{l>{2},  Fsf{2>); 

°l  normalization 
30  y  =  (yl  +  y2)  /  sqrt  (2)  ; 
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Listing  A. 15:  Code/DualTree/sfbDT.m 

1  function  y  =  sfbDT(lo,  hi,  sf) 


7,  Synthesis  filter  bank 


7 

USAGE : 

5  7 

y  = 

sf b (lo  ,  hi  ,  sf ) 

7 

INPUT  : 

7 

lo  - 

low  frqeuency  input 

7 

hi  - 

high  frequency  input 

7 

sf  - 

synthesis  filters 

10  7 

sf  (  : 

,  1)  -  lowpass  filter 

(even  length) 

7 

sf  (: 

,  2)  -  highpass  filter 

(even  length) 

7 

OUTPUT  : 

7 

y  - 

output  signal 

7,  See  also  afb 

15  7.  WAVELET  SOFTWARE  AT  POLYTECHNIC  UNIVERSITY,  BROOKLYN, 
7.  http  :  /  / taco  .  poly  .  edu/WaveletSoftware/ 


N  =  2* length ( lo )  ; 

L  =  length ( sf )  ; 

20  lo  =  upfirdn(lo,  sf ( : , 1 ) ,  2,  1); 

hi  =  upfirdn(hi,  sf ( :  ,  2)  ,  2,  1); 

y  =  lo  +  hi; 

y ( 1 : L -2 )  =  y(l:L-2)  +  y (N+ [1 : L-2] ) ; 
y  =  y  (  1  :  N)  ; 

25  y  =  cshift(y,  l-L/2); 


Listing  A.  16:  Code/DualTree/FSfarras.m 

1  function  [af ,  sf]  =  FSfarras 


7.  F arras  filters  organized  for  the  dual-tree 
7.  complex  DWT . 

5  7.  USAGE: 

7.  [af  ,  sf]  =  FSfarras 
7.  OUTPUT  : 
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•/.  af{i},  i  =  1,2  -  analysis  filters  for  tree  i 

1  sf{i},  i  =  1,2  -  synthesis  filters  for  tree  i 

10  1  See  farras,  dualtree,  dualfiltl. 

7.  WAVELET  SOFTWARE  AT  POLYTECHNIC  UNIVERSITY,  BROOKLYN,  NY 
7,  http  :  /  / taco  .  poly  .  edu/WaveletSoftware/ 


af  {  1 }  =  [ 


15 


20 


0 

-0 . 08838834764832 
0 . 08838834764832 
0 . 69587998903400 
0 . 69587998903400 
0 . 08838834764832 
-0 . 08838834764832 
0.01122679215254 
0.01122679215254 
0 


0 

-0 .01122679215254 
0.01122679215254 
0 . 08838834764832 
0 . 08838834764832 
-0 . 69587998903400 
0 . 69587998903400 
-0 . 08838834764832 
-0 . 08838834764832 
0 


25  ]  ; 


sf{l}  =  af { 1 } ( end : - 1 : 1 ,  :); 


af{2>  =  [ 

30  0.01122679215254 

0.01122679215254 
-0 . 08838834764832 
0 . 08838834764832 
0 . 69587998903400 
35  0.69587998903400 

0 . 08838834764832 
-0 . 08838834764832 
0 
0 


0 

0 

-0 . 08838834764832 
-0 . 08838834764832 
0 . 69587998903400 
-0 . 69587998903400 
0 . 08838834764832 
0 . 08838834764832 
0.01122679215254 
-0 .01122679215254 


40  ]  ; 


sf {2}  =  af {2} ( end : -1 : 1 ,  : ) ; 
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Listing  A. 17:  Code/DualTree/dualfiltl.m 


function  [af ,  sf]  =  dualfiltl 


7o  Kingsbury  Q-filters  for  the  dual-tree  complex  DWT 
7.  USAGE: 

%  [af ,  sf]  =  dualfiltl 

7.  OUTPUT  : 

7.  aflil,  i  =  1,2  -  analysis  filters  for  tree  i 

7.  sf{i},  i  =  1,2  -  synthesis  filters  for  tree  i 

7.  note:  af{2}  is  the  reverse  of  af{l} 

7.  REFERENCE: 

7«  N.  G.  Kingsbury,  "A  dual-tree  complex  wavelet 
7.  transform  with  improved  orthogonality  and  symmetry 
7.  properties",  Proceedings  of  the  IEEE  Int  .  Conf  .  on 
7.  Image  Proc  .  (ICIP),  2000 
7.  See  dualtree 

7.  WAVELET  SOFTWARE  AT  POLYTECHNIC  UNIVERSITY,  BROOKLYN, 
7«  http  :  /  / taco  .  poly  .  edu/WaveletSoftware/ 

7.  These  cofficients  are  rounded  to  8  decimal  places. 


af{l>  =  [ 

0 . 03516384000000 
0 

-0 . 08832942000000 
0 . 23389032000000 
0 . 76027237000000 
0 . 58751830000000 
0 

-0 . 11430184000000 
0 
0 


0 

0 

-0 . 11430184000000 
0 

0 . 58751830000000 
-0 . 76027237000000 
0 . 23389032000000 
0 . 08832942000000 
0 

-0 . 03516384000000 


1  ; 


af{2>  =  [ 


0  -0.03516384000000 
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45 


-0 . 11430184000000 
0 

0 . 58751830000000 
0 . 76027237000000 
0 . 23389032000000 
-0 . 08832942000000 
0 

0 . 03516384000000 


0 . 08832942000000 
0 . 23389032000000 
-0 . 76027237000000 
0 . 58751830000000 
0 

-0 . 11430184000000 
0 
0 


sf{l}  =  af { 1 } ( end : - 1 : 1 ,  :); 

sf {2}  =  af {2} ( end : -1 : 1 ,  : ) ; 
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